Aerodynamic Load Prediction on a Patrol Vessel Using Computational Fluid Dynamics

Abstract

Aerodynamic loads and moments on a naval patrol vessel are investigated using computational fluid dynamic simulations based on the OpenFOAM solver. After the initial turbulence, time, and grid dependency study, model scale simulations were performed for a wide range of inflow angles to predict aerodynamic forces and moments acting on the vessel at different heading conditions. For validation, model scale results were compared with wind tunnel data for similar hull forms. Finally, full-scale simulations were performed for a few cases to investigate possible scale effects on simulation results. The revealed scale effect turned out significant only for the yaw moment response. In this study, we aimed to produce reliable aerodynamic load data for the high-speed vessel, which is essential to developing reliable manoeuvring models. We conclude that Computational Fluid Dynamics is capable of providing reliable aerodynamic load predictions for high-speed vessels with sophisticated superstructures, in an economical manner.

1. Introduction

Although aerodynamic loads on surface displacement ships are typically viewed as secondary ones, their prediction is important for analysing the intact stability and simulating manoeuvring motion in wind [1,2]. Reliable information on aerodynamic loads can also be potentially used in the inverse problem of identification of the mathematical model parameters using full-scale ship data typically obtained in the presence of significant wind [3,4,5,6,7,8]. Another important application of the aerodynamic forces is to assess the ship speed loss due to wind and its effect on voyage duration and ship emissions [9,10,11].

As direct full-scale measurements of aerodynamic loads are practically impossible, these are typically estimated either through wind-tunnel tests with scaled models of the above-water part of the hull or using CFD codes. Several experimental or numerical results can be found in the available literature and part of these results are presented in the form of systematised databases suitable for a quick and rough estimation of aerodynamic forces and moments. Detailed wind tunnel test data for different vessel types were obtained by Owens and Palo [12], and somewhat later by Blendermann [13,14]. Andersen [15] provided some benchmark data for a post-Panamax container ship and made recommendations to make the container stack profiling smoother. As an alternative to using these databases of wind load data, approximate methods based on neural networks have also been adopted by Haddara and Guedes Soares [16].

With the development of computational fluid dynamics (CFD) and computational resources, several researchers have also started focusing on CFD for reliable prediction of wind loads. Hassan et al. [17] used CFD to investigate wind loads on a full-scale 2800 TEU container ship and discussed the optimization of stack arrangements. Koop et al. [18] investigated five different ship models and compared wind tunnel test results with those obtained with CFD. Luquet et al. [19] performed CFD simulations for an F70 frigate in heeled conditions and estimated the aero- and hydrodynamic coefficients. Wnȩk and Guedes Soares [20] conducted a CFD study of the wind loads on an LNG shuttle tanker in the vicinity of an offshore LNG production platform and compared the results with wind tunnel test data obtained by Wnȩk et al. [21]. Watanabe et al. [22] and Nguyen et al. [23] investigated wind resistance and moment acting on a 20,000 TEU container ship using CFD and discussed possibilities of reduction of the air drag. Majidian and Azarsina [24] investigated a 9000 TEU container ship following the Andersen [15] model. The same authors also investigated a 2748 TEU container vessel using CFD and proposed a wind resistance model based on statistical data [25]. Wang et al. [26] studied a cruise ship’s aerodynamics, experimentally and using CFD, and discussed design optimization possibilities through analysis of flow separation regions.

Apart from wind load studies, another group of researchers was involved in studies of ship air wakes, which are important for modelling helicopter landings and propagation of exhaust gases [27,28,29,30,31,32]. In particular, Linton and Thornber [33] presented a comprehensive review on CFD practices for ship air wake simulations, together with a discussion on a parameter dependency study and a method describing uncertainty quantification in associated turbulent flows. However, air wake studies are more focused on accurate descriptions and reproduction of the flow contours, including replication of the observed turbulence intensity and flow propagation. Whereas, in the case of load prediction, some compromise may be accepted as long as the average turbulence stress is well captured. Aerodynamic load studies are more focused on associated averaged Reynolds stress; thus, the total involved forces and moments, rather than the detailed flow field. As such, the study target and approach are quite different compared to air wake studies.

Although the studies performed for the determination of wind loads are numerous, they do not cover all existing types of ships and the available databases are not always sufficient for reliable predictions, which explains the necessity for individual studies contributing to the knowledge in ship aerodynamics. Practically, most of the existing aerodynamic load studies have been carried out in a model scale and the possible scale effect on forces and moments has rarely been discussed. The present study focuses on studying wind loads on a relatively small patrol ship of the Portuguese Navy with the help of the open-source Reynolds-averaged Navier–Stokes (RANS)-based CFD toolkit, OpenFOAM. While the main purpose of the study was to obtain a set of data on aerodynamic coefficients for systematically varying inflow angles for the model scale, the scale effect was also studied by means of a limited number of full-scale computations for a subset of the inflow angles. Similarly, the influence of the type of the involved turbulence model was investigated.

2. CFD Computational Layout

2.1. Ship Model

The ship under study belonged to the Portuguese Navy and is typically referred to as the Tejo class vessel. Historically, this design belonged to the family of Flyvefisken-class [34,35] patrol vessels primarily designed for the Royal Danish Navy.

The main particulars of the ship are length overall $L_{OA}=54\mathrm{m}$, beam $B=9\mathrm{m}$, and the full displacement draught $T=2.5\mathrm{m}$. Displacement of the vessel can vary from $320{\mathrm{m}}^3$ to $450{\mathrm{m}}^3$ and the maximum design speed reaches 30 knots. The image of the actual vessel (a) and the simplified model (b) used for simulations are shown in Figure 1. Since the purpose of the study was to estimate the involved aerodynamic loads and not the turbulence intensity, a simplified model was assumed sufficient for the purpose. Furthermore, being a multipurpose vessel, the deck arrangement of the vessel changed based on the intended function. Thus, the simplified model helped provide a generalised prediction for different deck arrangements or functionalities.

2.2. Mesh

For the study, single-phase simulations were performed only for the above-water part (freeboard) of the vessel with a somewhat simplified geometry of the superstructure. Thus, the mesh was generated only for the portion above the water, with the waterplane being a flat surface.

The simulation domain was generated following the general ITTC-2011 guidelines [36], with the inlet placed one ship length before the bow, the outlet placed two ship lengths after the stern, and the side boundaries at one ship length in each lateral direction from the sides of the hull (to the starboard and port). Although the ITTC guidelines are for hydrodynamic studies, a literature review by Linton and Thornber [33] revealed that the domain size has limited influence on aerodynamic simulations, unless the domain is too small. Furthermore, flow contours from the simulations did not reveal any obvious boundary interactions. Thus, the adopted domain size was considered sufficient.

The bottom of the domain was modelled as a rigid wall fixed at the waterplane, and the upper boundary (atmosphere) was located at about half the ship length above the bottom boundary. This upper boundary and the side boundaries were treated as local symmetry planes. The symmetry boundary condition was applied at a mirror surface, where the fluxes and normal components of all variables across the symmetry were set to zero. The symmetry in this case works as a non-reflecting boundary condition that ignores the interaction of the flow field with the sides and the top. The pressure boundary condition at the inlet is a zero-gradient (the OpenFOAM function zeroGradient: Neumann boundary condition), and the outlet boundary supplies a constant pressure condition (fixedValue: Dirichlet boundary condition). Regarding the velocity field, the inlet was set as a constant velocity condition (fixedValue), and the outlet was set to be a generic outflow condition with zero return inflow (inletOutlet). Initial values for the turbulence parameters were calculated following the Reynolds number [37]. The OpenFOAM built in a rough wall function was used to reduce the mesh dependency near the wall to capture the boundary layer. The airflow at the inlet was assumed to be uniform. The domain size and the domain boundaries are shown in Figure 2.

The simulation domain was created using the blockMeshDict utility from OpenFOAM, which generates a block of a defined dimension with structured grids. The block was refined consecutively four times to obtain a higher mesh resolution around the vessel. Separate refinement blocks were used for the hull up to the deck and for the superstructure above. This was done to avoid unnecessary refinement in empty areas and to reduce the total mesh size. Finally, the hull was integrated into the domain using the snappyHexMesh utility, which uses the Boolean approach for merging the solid object into the grid space. Five layers around the hull were generated using the function snappyHexMesh to ensure better capturing of the viscous resistance. The refinement area and the hull position were adjusted by changing the inflow angle while keeping the simulation domain unchanged. In the domain, the ship’s longitudinal direction was in a positive $x$-axis (stern to bow), the lateral section was in the $y$-axis and the vertical direction was in the positive $z$-axis.

For the model scale simulations (1:10 scale), an average mesh resolution of 3.95 million cells was used, with a non-dimensional wall distance, $y^{+}\approx 87$, where $y^{+}$ was used to describe the height of the first grid element next to the wall. Further details about the model scale mesh resolution are provided in the verification study section. The full-scale simulations were performed with a similar mesh topology with a substantially higher mesh resolution of around 80 million cells. A general image of the simulation mesh and refinement around the vessel is shown in Figure 3.

2.3. The Solver and Computational Resource

As previously mentioned, the study is based on the OpenFOAM library, whose solvers have been elaborately described by Weller et al. [38] and Jasak [39]. For the presented study, the OpenFOAM version 2006, managed by the ESI group, (https://www.openfoam.com/openfoam-is-open-source) (accessed on 15 August 2021) was used.

The governing equations for the solver applied here are the Reynolds-averaged Navier–Stokes (RANS) equations coupled with the continuity equation. In the vector–tensor form for an incompressible, turbulent, single-phase unsteady flow, these equations are written as:

$$\begin{gathered} \frac{\partial v}{\partial t}+v·\nabla v-\nabla · \left(\nu \nabla v\right)=\frac{1}{\rho }\nabla p, \\ \nabla ·v=0 \end{gathered}$$

(1)

Here, $t$ is the time, $\nabla$ is the Hamilton operator describing the gradient or divergence,$v$ is the fluid velocity, $\rho$ is the fluid density, $p$ is the pressure, $\nu$ is the effective viscosity. The finite volume method (FVM) is used to discretise the governing equations. The pressure–velocity coupling is handled through the PIMPLE algorithm (a combination of PISO and SIMPLE algorithms). OpenFOAM incorporates several turbulence models, and for the present paper, the SST $k-\omega$ model is used for most computations. All computations are in a single phase and the free surface was modelled as a flat surface with no friction.

For the simulations, a desktop computer with Intel i9 processors with 18 cores, 128 GB RAM, and a 2TB SSD disk was used. The model scale simulations were performed using 12 processors in parallel, while the full-scale simulations were performed using 18 processors. Each model scale simulation was run for at least 25 simulation seconds (30 flow-through times) to reach stable outputs; that took roughly 40 physical hours. As for the full-scale simulations, each case was simulated for roughly 40 s of simulation time (15 flow-through times); that took roughly 650 h. For this case, the simulation time was reduced due to large resource dependency and because the solution converged (steady response) earlier due to smaller time steps.

3. Results

The OpenFOAM PIMPLE (single-phase) solver was used for both the model- and full-scale simulations. Only the slightly simplified above-water vessel geometry was involved and it was assumed symmetric with respect to the centre plane. Simulations were performed on the vessel model remaining fixed at each value of the inflow angle, assuming zero heel and fixed trim. The air at the inlet was assumed to have constant velocity without any pre-defined gust or flow oscillation.

A 1:10 model (5.4 m in length) was used for the model scale simulations. Simulations were performed at the inflow velocity of 20 m/s for the full-scale and 6.3 m/s for the model scale. The corresponding values of the Reynolds number based on the overall ship and model length were $7.12\times {10}^7$ and $2.28\times {10}^6$, respectively, assuming an air temperature of 20 °C.

In experimental aerodynamic studies, the focus is mostly given to maintaining Reynolds and Strouhal similarity, since it is easier and less time-consuming to reach high Reynolds numbers in laboratories [38]. However, a compromise had to be made in the CFD study considering resource limitations. Even though boundary layer and, thus, frictional resistance, have limited impacts on the total aerodynamic load, maintaining a reasonable y⁺ value near the wall also ensures that the meshing follows reliable guidelines and mesh resolution around the hull is sufficient to properly capture the vortices and pressure forces. To maintain Reynold’s similarity, the mesh requirement for the model scale becomes even higher compared to a full-scale simulation. Thus, to maintain the practicality of the study, Froude similarity has been applied, ignoring the wind tunnel practice of maintaining Reynold’s similarity.

The CFD computations were performed in the following steps:

1.

Turbulence model dependency study for a mesh with 3.95 million cells using three different approaches, by modelling with SST K-Omega and K-Epsilon two-equation turbulence models, and by resolving with the Delayed Detached Eddy Simulation.

2.

Time convergence study for three inflow angles (0°, 45°, 90°, 135°, and 180°) using a mesh with 3.95 Mcells and the selected (SST K-Omega) turbulence model.

3.

Grid and time step dependency study for a 45° inflow angle using three mesh resolutions and the selected turbulence model.

4.

Study of aerodynamic forces and moments encountered by the vessel while facing wind from different inflow angles. In total, the resulting database includes wind loads for 24 inflow angles ranging from head to stern flow.

5.

Study of scale effects on aerodynamic load prediction by simulating five selected cases in full-scale and comparing with model scale results.

The results are presented in the form of several visualisations (flow contours and pressure distributions) to show some peculiarities of the flow, by plotted time histories reflecting convergence studies and by plotted and tabulated results for the aerodynamic forces, moments, and their normalised values. All time histories are shown from 2 to 20 s, to show the initial oscillation before reaching a steady response. A longer time history was avoided to make the oscillations clearly visible, which becomes difficult to represent with longer time histories (especially since some figures had to be presented in smaller dimensions). Primary results for the forces were obtained in the velocity axes in form of the drag $F_x$ and lateral force $F_y$, and were used directly in the convergence studies. The final study results are represented in body axes by means of evident transformations.

Direct validation of the results for the ship under study has not been possible but comparisons with somewhat similar shapes from Blendermann’s database [13] were performed.

3.1. Preliminary Studies

3.1.1. Turbulence Models

The interaction between the vessel’s upper structures and the high-velocity flow is a highly turbulent phenomenon. For accurate modelling of the forces involved in the interaction, reliable modelling of the turbulence is essential. In air wake studies, mostly the Large Eddy Simulation (LES) or Detached Eddy Simulation (DES) are used [33]. However, LES and DES are computationally too expensive and normally not practical for series computations. In the present study, two simpler and commonly known turbulence models and a numerical approach are tested. These are the two-equation turbulence models (SST K-Omega [40] and K-Epsilon [41,42]) and the Delayed Detached Eddy Simulation (DDES) [43]. Among the three approaches, the turbulence models model the turbulence, whereas, the DDES attempts to resolve the turbulence.

Simulations were run for the 45°, 90°, and 135° air drift angles with the three approaches and with the same mesh resolution of 3.95 million. All simulations were performed for at least 30 flow-through times.

The body axis-based non-dimensional surge and sway force coefficients $C_{XA}$ and $C_{YA}$ (as defined in Section 3.3) and roll and yaw moment coefficients $C_{KA}$ and $C_{NA}$ (as defined in Section 3.3) are presented in

Table 1. The surge and sway force and heel and yaw moment results for the aerodynamic load-predicted using model scale simulations in OpenFOAM.

Drift Angle (Deg)	Turbulence Model	${\mathit{C}}_{\mathit{XA}}$	${\mathit{C}}_{\mathit{YA}}$	${\mathit{C}}_{\mathit{KA}}$	${\mathit{C}}_{\mathit{NA}}$
45	SST K-Omega	−0.22	0.85	−0.522	0.16
	K-Epsilon	−0.19	0.85	−0.515	0.16
	DDES	−0.13	0.92	−0.564	0.17
90	SST K-Omega	−0.04	1.10	−0.599	0.06
	K-Epsilon	0.03	1.14	−0.587	0.06
	DDES	−0.12	1.17	−0.627	0.06
135	SST K-Omega	0.36	0.93	−0.458	−0.06
	K-Epsilon	0.29	1.00	−0.484	−0.08
	DDES	0.29	0.98	−0.471	−0.06

. The results show somewhat oscillating behaviours with different models predicting higher values in different cases. Nevertheless, the overall differences among the results remain minimum for most cases. Because of the absence of experimental data, it is difficult to conclude which model predicts better. However, the relatively low mesh resolution for the DDES cases might have contributed to the observed over-predictions.

Considering that the same mesh resolution was applied, the time required for performing the simulations was roughly similar. However, in general, simulations with the K-Epsilon model were faster compared to the other two, and the model with a relatively slower solution was SST K-Omega.

The force and moment time history results from the simulations are presented in Figure 4 and Figure 5, respectively. The time histories shown are based on the velocity axis. The figure shows drag $F_x$ and transverse force $F_y$ histories, and the roll, pitch, and yaw moments ($M_x=K$, $M_y=M$, and $M_z=N$) time histories.

The images show that for all three turbulence models, maximum oscillation in the force and moment was observed in the 90° drift simulation. For all drift angles, the K-Epsilon model showed the most stable results with almost a linear flat line for forces and moments. This might suggest that the K-Epsilon turbulence model oversimplifies the turbulence in this case and provides average stress results instead of considering the oscillating turbulent phenomenon. The SST K-Omega model shows better capturing of turbulence, which is represented by the repeated oscillations in the time history within certain bounds. As for the DDES model, it shows the highest oscillation for the force and moment recordings compared to the other two models. In theory, DDES is supposed to better represent the turbulence phenomenon, which explains the higher oscillation. However, DDES requires a substantially high mesh resolution, both near and far from the body (vessel), and low CFL conditions to properly capture the turbulent behaviour of the flow. Considering the applied mesh resolution, DDES results would most probably reveal high uncertainty in results, thus reducing the reliability of the study. Moreover, the application of a higher mesh resolution is not feasible considering the large number of simulations planned for the study.

Turbulence modelling is a complex issue and assessing different turbulence models without direct comparison with experimental data is challenging. Different models might be better suited for different types of fluid–structure interactions and different regions. As such, making conclusive statements about the better suitability of different turbulence models for a particular study is very difficult with an extensive investigation. As mentioned earlier, the purpose of the present study was not to observe or measure the turbulence intensity involved in the interaction, but rather to study the forces and moments involved. As such, the SST K-Omega turbulence model was selected for the study, which seemed most reliable for this study after the investigation.

3.1.2. Time Domain Settling

CFD simulations in the time domain are performed until some steady values of the forces are obtained or a regime of somewhat stationary oscillations is reached. The latter may occur when non-streamlined shapes are placed into a steady flow. In those cases, averaging over the last 10 s (12 flow-through times) of simulation time, after reaching a repeated steady response, is performed to obtain the mean values of the loads. For full-scale simulations, a larger flow-through time is considered for averaging. All simulations after the initial turbulence study are performed using the SST K-Omega turbulence model.

The drag and the transverse force time histories for five simulation cases (0°, 45°, 90°, 135°, and 180°) are shown in Figure 6. The roll, pitch, and yaw moment histories are shown in Figure 7, for the same cases. However, as before, the time histories shown are based on the velocity axis.

Although the time history for the first 20 s of the simulation (24 flow-through times) is shown in the figures, simulations were run for a longer period, especially for the cases with high oscillation (e.g., the 90 deg case). The figures show that relatively high oscillation is observed for the beam wind case. The oscillation can be partly explained by the high turbulence associated with this case, as can be seen in Figure 8.

3.1.3. Grid Convergence Study

Next, a systematic verification study was performed for the simulation settings and mesh configuration to assess the associated numerical uncertainty. Three different mesh configurations were used to perform a systematic verification study following the processes suggested by ITTC-2008 [44] guidelines, namely the factor of safety or the $F_s$-based approach [45]. However, contrary to the ITTC guidelines, instead of conducting an independent time and grid-based study, a constant Courant–Friedrichs–Lewy (CFL)-based discretization approach is used with equal refinement for both the grid and the time step. The justification for this approach has been discussed by Islam and Guedes Soares [46,47] and Wang et al. [48]. The details of the mesh resolutions used for the model scale uncertainty study are shown in Table 1. The table shows the total cell number in the domain and the general cell dimensions for each mesh in $x,y$, and $z$ directions. The minimum layer thickness represents the minimum cell dimension in the direction normal from the hull surface). The layer is used to ensure the minimum required cell size near the hull surface to meet the $y^{+}$-criterion for boundaries with wall function (30 < $y^{+}$ < 300, following the law of wall). The refinement ratio is defined as: $r_{n+1,n}=h_{n+1}/h_n$, where $h$ is the dimension of a cell in any particular direction ($x,y$, and $z$), and the subscript $n=1,2$ is the mesh number, according to

Table 2. Mesh resolutions used for the verification study for the model scale simulations.

Mesh	Total Number of Cells	Dimensions x (m)	Of y (m)	Cells $\mathit{z}$	Minimum Layer Thickness	Non-Dimensional Wall Distance, y+	Coarsening Ratio
1	$8.7\times {10}^6$	0.01953	0.01875	0.01875	0.00375	65	1.00
2	$3.95\times {10}^6$	0.02500	0.02500	0.02344	0.00500	87	1.30
3	$2.25\times {10}^6$	0.03125	0.03125	0.03125	0.00625	108	1.25

.

For the uncertainty study, all computations were performed for the inflow angle of 45 degrees. The uncertainty was quantified for the surge ($X_A$) and sway ($Y_A$) force, and roll ($K_A$) and yaw ($N_A$) moments. The results are shown in

Table 3. Uncertainty analysis for the aerodynamic forces and moments.

Property		${\mathit{X}}_{\mathit{A}} \left(\mathbf{N}\right)$	${\mathit{Y}}_{\mathit{A}} \left(\mathbf{N}\right)$	${\mathit{K}}_{\mathit{A}} \left(\mathbf{Nm}\right)$	${\mathit{N}}_{\mathit{A}} \left(\mathbf{Nm}\right)$
Simulation results	$S_1$ (fine)	−3.73	50.18	−16.45	53.25
	$S_2$ (mid)	−3.65	49.95	−16.68	50.20
	$S_3$ (coarse)	−4.63	49.84	−16.47	48.48
Refinement ratio	$r_{21}=h_2/h_1$	1.3	1.3	1.3	1.3
	r32 = h3/h2	1.25	1.25	1.25	1.25
Difference of estimation	${\epsilon }_{21}=S_2-S_1$	0.0778	−0.233	−0.233	−3.050
	${\epsilon }_{32}=S_3-S_2$	−0.985	−0.114	0.208	−1.720
Convergence	${\epsilon }_{21}/{\epsilon }_{32}$	−0.079	2.050	−1.122	1.773
Order of accuracy	$p$	11.25	1.96	0.41	1.45
Extrapolated values	$S_{ext}^{21}$	−3.73	50.53	−14.39	59.84
	$S_{ext}^{32}$	−3.56	50.16	−18.85	54.70
Approximate relative error	$e_a^{21}$	−0.0209	−0.0046	0.0142	−0.0573
	$e_a^{32}$	0.270	−0.0023	−0.0125	−0.0343
Extrapolated relative error	$e_{ext}^{21}$	−0.00115	−0.0069	0.1428	−0.1101
	$e_{ext}^{32}$	0.0245	−0.0041	−0.1151	−0.0823
Grid convergence index (GCI)	GCI21	−0.0014	−0.0086	0.15616	−0.15467
	GCI32	0.0298	−0.0052	−0.1626	−0.1121
Corrected uncertainty	$U_{c1}$	0.0288%	0.1729%	3.1231%	3.0933%
	$U_{c2}$	0.5968%	0.1039%	3.2522%	2.2421%

, where $S_{1,2,3}$ stands for the computation results, per se. The absolute value of the convergence factor indicates whether the solution is converging (<1) or diverging (>1). For predicting the apparent order of accuracy, $p$, the following equation was used by Celik et al. [45]:

$$\begin{array}{l}p=\frac{1}{\ln (r_{21})}\left|\ln \left|\frac{{\epsilon }_{32}}{{\epsilon }_{21}}\right|+q(p)\right|,\\ q(p)=\ln \left(\frac{r_{21}^p-s}{r_{32}^p-s}\right),\\ s=\mathrm{sgn}\left(\frac{{\epsilon }_{32}}{{\epsilon }_{21}}\right).\end{array}$$

(2)

Celik’s paper also proposed the following extrapolated values for the solutions:

$$s_{ext}^{21}=\left|\frac{r_{21}^p{S}_1-S_2}{r_{21}^p-1}\right|$$

(3)

and the estimates of the errors and the extrapolated errors are:

$$\begin{array}{l}{e}_a^{21}=\left|S_1-\frac{S_2}{S_1}\right|,\\ e_{ext}^{21}=\left|S_{ext}^{21}-\frac{S_1}{S_{ext}^{21}}\right|\end{array}$$

(4)

The Grid Convergence Index (GCI), which is the uncertainty estimate $U_i$, is defined as a product of the Richardson normalised discretization error (${\delta }_{RE,1}^{*}$) by $F_s$ (factor of safety). The value $F_s$ is suggested to be 1.25 for the systematic parameter refinement study with at least three inputs, or 3 for simple convergence studies with two values of the input parameter.

$$GC{I}_{21}=\frac{1.25\hspace{0.17em}{\delta }_{21}^a}{r_{21}^p-1}=F_s\left|{\delta }_{RE,1}^{*}\right|$$

(5)

$${\delta }_{21}^a=\left|\frac{S_2-S_1}{S_1}\right|$$

(6)

Following the ITTC-2008 guidelines, the corrected uncertainty ($U_c$), is defined as:

$$U_c=(F_s-1)\left|{\delta }_{RE,1}^{*}\right|$$

(7)

The uncertainty study shows (Table 3) oscillatory convergence for the surge force and the roll moment and divergence for the sway force and yaw moment results. In general, for high Reynolds number simulations with complex structures, the involved turbulence makes it very difficult to have monotonously converging results. Furthermore, the study angle of 45 degrees represents a very turbulent condition. Nevertheless, in all the cases, the estimated corrected uncertainty is less than 3%, except for the rolling moment. Direct validation of the results could not be attempted due to the absence of experimental data for the same ship model. Nevertheless, a comparison is shown in Section 3.4 for the simulated results and experimental results for similar hull forms to assess the reliability of the predictions.

3.2. Flow Field Visualization Results

A detailed study of the aerodynamic loads on the Tagus class vessel was performed using a large number of model scale simulations. Single-phase simulations were performed for the vessel facing wind at 6.3 m/s velocities from different directions. Assuming symmetry of the above-water part of the vessel, the wind inflow direction was changed from the ship bow to the stern with the 10 deg interval, but some additional values of the inflow angle were also investigated.

Initially, flow field visualizations from the turbulence model study were presented in Figure 8, for 45- and 90-degree drift angles. The flow field contours further confirmed the findings discussed in Section 3.1.1, based on the force and moment results. The figure shows higher chaotic behaviour of the flow in the case of DDES and relatively more streamlined flow for the K-Epsilon model.

Next, flow field simulation results are presented for some specific airflow angles, performed using the defined mesh resolution (3.95 mil) and turbulence model (SST K-Omega), to understand how wind interacts with the vessel above the water structure and how forces and moments are influenced.

A flow field visualization for the wind interacting with the model scale vessel for the head and the following wind is shown in Figure 9. The figure shows a vortex formation after (for headwind) and before (for the following wind) the cabin structure. The streamlines show changes when the inflow angle is changing, and the colour on the hull (above the water part) shows pressure distribution due to wind interaction. The high-pressure region at the front of the superstructure in the head waves indicates relatively higher wind resistance.

To further illustrate the results, Figure 10 shows the flow stream interacting with the vessel at different angles. The contours show the change in wind direction and vorticity formation after interacting with the vessel. The streamlines represent the change in free stream velocity and pressure distribution on the hull superstructure (as shown). The turbulent airflow in the case of the beam sea indicates a higher roll motion. The pressure distribution, in this case, is also different because of local recirculation regions. Notable turbulence is also observed in 45° and 135° cases.

After model scale simulations, full-scale simulations were performed from the head to the following wind with a 45° interval. The number of simulation cases was reduced considering the substantial amount of computational resources required for each full-scale simulation.

Pressure distribution on the hull form (above water) and free stream flow around the hull are shown in Figure 11 for the full-scale simulations. The figures show notably higher turbulence compared to the model scale simulations. The variation mostly comes from the difference in the Reynolds scale, since in full-scale simulations, turbulence is fully developed. This enhanced turbulence would lead to a higher roll moment prediction in a full-scale simulation compared to the model scale.

3.3. Force and Moment Results

As mentioned earlier, the forces and moments were brought to the right-handed Cartesian body frame with the origin in the intersection of the midship plane, the centre plane, and the waterplane; with the $x$-axis pointing forward, $y$-axis—to the starboard, and the $z$-axis—downwards. The aerodynamic forces of surge $X_A$ and sway $Y_A$, and moments of roll $K_A$ and yaw $N_A$ defined in body axes can be represented as:

$$\begin{array}{l}{X}_A=C_{XA}\left({\beta }_A\right)\frac{\rho V_A^2}{2}{A}_T,\hspace{1em}{Y}_A=C_{YA}\left({\beta }_A\right)\frac{\rho V_A^2}{2}{A}_L,\hspace{1em}\\ K_A=C_{KA}\left({\beta }_A\right)\frac{\rho V_A^2}{2}{A}_L{L}_{OA},\hspace{1em}{N}_A=C_{NA}\left({\beta }_A\right)\frac{\rho V_A^2}{2}{A}_L{L}_{OA},\end{array}$$

(8)

where $C_{XA}, C_{YA}, C_{KA}, C_{NA}$ are the corresponding aerodynamic force/moment coefficients; ${\beta }_A$ is the air drift or inflow angle,$A_T$ and $A_L$ are, respectively, the frontal (transverse) and the lateral projected area of the above-water hull; $L_{OA}$ is the length overall.

The model scale simulation results of aerodynamic loads are shown in

Table 4. The computed aerodynamic forces/moments in the model scale and their coefficients.

Drift Angle	${\mathit{X}}_{\mathit{A}} \left(\mathbf{N}\right)$	${\mathit{Y}}_{\mathit{A}} \left(\mathbf{N}\right)$	${\mathit{K}}_{\mathit{A}} \left(\mathbf{Nm}\right)$	${\mathit{N}}_{\mathit{A}} \left(\mathbf{Nm}\right)$	${\mathit{C}}_{\mathit{X}\mathit{A}}$	${\mathit{C}}_{\mathit{Y}\mathit{A}}$	${\mathit{C}}_{\mathit{K}\mathit{A}}$	${\mathit{C}}_{\mathit{N}\mathit{A}}$
0	−7.90	0.17	−0.07	1.00	−0.48	0.00	−0.002	0.00
10	−9.19	9.33	−2.39	15.94	−0.56	0.16	−0.075	0.05
20	−8.68	24.02	−7.49	29.93	−0.53	0.41	−0.234	0.09
30	−6.75	39.57	−12.41	43.33	−0.41	0.68	−0.388	0.14
40	−4.71	45.88	−15.05	51.10	−0.29	0.78	−0.471	0.16
45	−3.65	49.95	−16.70	50.20	−0.22	0.85	−0.522	0.16
50	−2.82	52.97	−17.44	49.76	−0.17	0.90	−0.545	0.16
60	−3.71	56.42	−16.75	43.31	−0.23	0.96	−0.524	0.14
70	−3.38	62.12	−17.53	37.22	−0.21	1.06	−0.548	0.12
80	−3.15	66.06	−19.85	27.79	−0.19	1.13	−0.621	0.09
85	0.22	63.21	−19.62	20.00	0.01	1.08	−0.614	0.06
90	−0.62	64.51	−19.15	18.41	−0.04	1.10	−0.599	0.06
95	−2.40	66.04	−19.21	7.71	−0.15	1.13	−0.601	0.02
100	−2.15	64.86	−20.09	4.56	−0.13	1.11	−0.628	0.01
110	0.63	63.21	−17.03	−9.72	0.04	1.08	−0.533	−0.03
120	2.19	62.15	−16.69	−10.20	0.13	1.06	−0.522	−0.03
130	5.00	55.55	−15.20	−15.32	0.30	0.95	−0.475	−0.05
135	5.89	54.58	−14.66	−18.81	0.36	0.93	−0.458	−0.06
140	6.37	51.20	−13.93	−21.91	0.39	0.87	−0.436	−0.07
145	7.70	46.99	−12.60	−22.67	0.47	0.80	−0.394	−0.07
150	7.73	41.70	−11.20	−21.95	0.47	0.71	−0.350	−0.07
160	8.83	28.57	−6.99	−14.12	0.54	0.49	−0.219	−0.04
170	9.97	13.03	−2.80	−7.38	0.61	0.22	−0.087	−0.02
180	8.25	0.20	−0.02	−0.27	0.50	0.00	−0.001	0.00

. The results show both the dimensional values for the model scale simulations and the non-dimensional coefficient.

The negative values represent the drag. According to the results, the vessel faces a maximum air drag at bow wind, which keeps decreasing as it approaches the beam wind conditions. After the beam, the drag continues for a while, and then it starts facing a lift force from the following wind. The sway force coefficient shows a gradual rise from the headwind to beam condition, reaches the maximum at the beam, and then starts decreasing again as it approaches following wind conditions. The roll moment results show that the maximum moment is observed close to the beam wind condition. As for the moment, the result shows notable asymmetry around the beam condition. The yaw moment observed is significantly higher for headwind conditions compared to the following wind cases.

The full-scale simulation results are shown in

Table 5. The computed aerodynamic forces/moments in full-scale and their coefficients.

Drift Angle	${\mathit{X}}_{\mathit{A}} \left(\mathbf{N}\right)$	${\mathit{Y}}_{\mathit{A}} \left(\mathbf{N}\right)$	${\mathit{K}}_{\mathit{A}} \left(\mathbf{Nm}\right)$	NA (N-m)	${\mathit{C}}_{\mathit{X}\mathit{A}}$	${\mathit{C}}_{\mathit{Y}\mathit{A}}$	${\mathit{C}}_{\mathit{K}\mathit{A}}$	${\mathit{C}}_{\mathit{N}\mathit{A}}$
0	−7396.12	−174.56	700.00	−1227	−0.44	0.00	0.00	0.00
45	−4201.27	48,412.00	−210,030.51	535,125	−0.25	0.82	−0.65	0.17
90	−1338.89	64,698.98	−251,004.00	135,944	−0.08	1.10	−0.78	0.04
135	5146.32	57,890.83	−221,197.14	−264,855	0.31	0.98	−0.69	−0.08
180	9579.00	−272.46	1490.47	−2939.8	0.58	0.00	0.00	0.00

with the predicted forces, moment data, and coefficients. The dimensional values from the full-scale simulations can serve as a reference to understand the encountered aerodynamic loads by the vessel.

The non-dimensional results for forces and moments from both the model scale and full-scale simulations are also shown in Figure 12 and Figure 13.

Relative comparison between the full-scale and model scale results (Figure 12 and Figure 13) suggests that the non-dimensional force results agree quite well with each other. Although slight over-predictions are observed in the case of model scale results for surge force, sway force results agree very well. Relatively notable variations are observed in the case of moment prediction, especially at the 135° drift angle.

The results suggest that the model scale simulations slightly underpredict the yaw moments in most cases. However, in the case of rolling moments, a notable difference is observed. This might be attributed to the fully developed turbulence in the actual Reynolds number, which was missing in the model scale. The variation, in this case, is mostly related to the difference in the Reynolds number (scale effect).

However, some contributions in the deviation also come from the simulation meshes, since the $y^{+}$ value for the model scale and full-scale simulations could not be maintained at the same level. Nevertheless, the overall difference between the model and full-scale results remains limited. This suggests that aerodynamic load coefficients show a minor impact from the simulation scale factor. This is mostly because aerodynamic loads are pressure dominant with a minor contribution from viscous forces. As such, scaling of the results mostly depends on Reynolds scaling and does not require empirical formulation for scaling of viscous forces, which often provide over prediction [49].

3.4. Comparison of Results with Blendermann Models

Direct validation of the results was not possible due to the absence of experimental or sea-trial data. Nevertheless, a general comparison with the results for a Destroyer (DES0101BN), a research vessel (RES0101BN), and a speedboat (SPE0102BN) model from Blendermann [13] wind load dataset suggests that, in general, OpenFOAM overpredicts the sway force and yaw moment.

The vessel outlines are shown in Figure 14. The navy destroyer (DES0101BN) is a 133.8 m vessel, with a breadth of 13.40 m, and with an above-water lateral area of 1342.07 m². The experiment was performed on a 1:100 scale. The research vessel had an overall length of 82.05 m, a breadth of 13.50 m, and the reference area was 670.27 m². The speedboat (SPE0102BN) had an overall length of 53.60 m, a breadth of 9.20 m, and a reference area of 317.63 m². As can be observed, the parameters and free-board geometry of the vessels considered here are notably different from the Tejo class vessel. Nevertheless, it provides a general understanding of the results.

The comparisons are shown in Figure 15 and Figure 16. Figure 15 shows a comparison for surge and sway forces, whereas Figure 15 shows a comparison for roll and yaw moments. The surge force coefficient in Figure 13 was nondimensionalised using the longitudinal windage area instead of the transverse area used in Figure 12. The differences in the responses seem to be quite adequate regarding notable differences in the shape of the above-water part of the hull.

4. Conclusions

In the study, RANS simulations for a navy vessel were performed to predict the aerodynamic drag of the vessel at different heading angles. Both model scale and full-scale simulations were performed using OpenFOAM, and drag coefficients were predicted, ignoring the free surface effect. Initially, a turbulence model dependency and a verification study were performed for model scale simulations to understand related dependencies. Next, 24 model scale simulations were performed for various heading angles. Finally, full-scale simulations were performed for five heading positions. Direct validation of the results was not possible due to the absence of experimental or sea trial data. However, a relative comparison of results with similar hull forms confirmed the reliability of the results.

The study provides surge and sway force and roll and yaw moment coefficients for the aerodynamic loads on Tejo class vessels. The result curves show reasonable agreement with similar numerical studies available in the literature. The model scale and full-scale simulations show acceptable agreement in most cases. The comparison of the model and full-scale results show that the scale effect is insignificant for the surge and sway forces and more significant for the yaw moment. Reasons for such differences remain subject to further investigation.

In general, following the dependency studies and comparison with available data, it may be concluded that CFD is well capable of providing reliable data for the aerodynamic loads acting on the studied vessel, which can be used in further studies to improve the manoeuvring trajectories of the vessel under the action of wind. The paper contributes by providing reliable quantitative predictions of aerodynamic loads and moments for a high-speed vessel. It also demonstrates where scale effects are more prominent. However, the paper adopts a relatively low mesh resolution and skips detailed capturing of turbulence intensity and propagation, which are important for several operational aspects of a navy ship. Thus, to assess the difference, a more detailed approach following the air wake study group for simulators might be a future study target.