Numerical simulation of container ship in oblique winds to develop a wind resistance model based on statistical data

ABSTRACT The forces that wind exerts on a merchant ship, especially for a container ship that carries high stacks of cargo on board, can be very large; therefore, more accurate simulations than presently available in the literature on the aerodynamics of a container ship are necessary for a safe and efficient shipping. This research tries to demonstrate a newer level of accuracy in capturing the aerodynamic properties of a typical full-stack container ship. Simulations setup was verified against US National Advisory Committee for Aeronautics (NACA) test data for the Akron airship. These are unique data that can verify the accuracy of the Computational Fluid Dynamics (CFD) prediction for pressure distributions over a three-dimensional (3D) object in nonzero wind angles. The numerical simulations of the container ship, which have been performed by ANSYS-CFX™, focus on wind angles of attack (AOAs) up to 40°, obtained from a voyage data record of a container ship throughout a rather long voyage. A regression equation for air resistance is derived within this range, and a set of voyage data for wind speed and directions is also analyzed to check the applicability of the results. Innovative features of air–sea interface modeling have been included in this report.


Introduction
The new generation of large container ships, such as post-Panamax size, can load eight tiers of containers in height, more than 30 rows in length, and up to 23 bays across the breadth on deck ("FleetMon," AIS live vessel tracking community 2018). A fully loaded container ship can load more than 50% of containers above its deck, and this differentiates it from other kinds of merchant vessels. Consequently, projection areas of such ships above the water surface are tremendous. In a container ship with a large windage area, wind resistance comprises a significant portion of the total resistance, up to 10% (Minsaas and Steen 2008).
In the past few decades, many studies have been carried out to investigate the effect of wind load on a container ship with large windage area from different aspects. Studies by Andersson (1978) and van Berlekom (1981) are two earlier reports. According to Andersson (1978), the induced resistance from the increased rudder angle plays an insignificant role, while according to van Berlekom (1981), it can be as large as the longitudinal force for stronger winds. Generally, the longitudinal wind force is of the greatest importance for predicting the propulsion force. Its share in the total resistance has been discussed by van Berlekom (1981) and Aage (1968), who stated that wind resistance rarely contributes more than 10% of the total resistance. For the first time, it was noted by Blendermann (1997) that the wind forces are influenced by container stacking configuration. Blendermann stated that uneven bay heights increase the wind resistance. A ship with randomly stacked containers on deck compared with a fully stacked ship experiences significantly higher longitudinal force, smaller transverse force, smaller rolling moment, and smaller yaw moment (Blendermann 1997). In Andersson (1978), the individual configurations have been more thoroughly described and discussedthe containers on the fore deck play a more important role in wind resistance than the aft deck, and random container configuration can significantly increase the longitudinal force. Finally, it was concluded that the configuration with full load on the fore deck and streamlining of the aft deck is the most favorable among the tested configurations. Andersen (2012aAndersen ( , 2012b investigated the wind forces acting on a 9,000+ Twenty Foot Equivalent Unit (TEU) post-Panamax container ship using a 1:450 scale model in an extensive series of container configurations on deck through wind tunnel tests. The results stated that a minimum drag was experienced when the containers were loaded on the aft and fore decks in streamlined configuration. Anderson performed the experimental wind tunnel model tests for different stacking configurations while longitudinal force and yaw moment were investigated.
The second study supported the previous results along with more details. Hassan, White, and Ciortan (2012), using the computational fluid dynamics approach, investigated the effect of container configuration on the wind resistance of a 2800 TEU container ship 221.65 m in length. They reported the longitudinal and lateral forces along with pressure contours and velocity streamlines. Fourteen different stacking configurations were investigated for wind directions from 0°t o 90°. Finally, the alteration of wind resistance and its effect on fuel consumption and gas emissions were investigated. Watanabe et al. (2016) carried out CFD simulations and measurements of wind resistance and moment acting on a scale model of a 20,000 TEU container ship. Nguyen et al. (2016) studied the longitudinal and transverse forces and yaw moment acting on the same scaled model using the CFD method. On the basis of numerical results, the authors proposed several methods to reduce the air resistance in head wind by using a winglet, reduction of the front containers, bow covers, and back guides. Wnek et al. (2015) also performed experiments on models of a liquefied natural gas carrier (LNGC) and a floating platform. They assessed the effect of wind loads and the induced motions at several different arrangements of the two vessels with respect to each other. Janssen, Blocken, and van Wijhe (2017) reviewed Andersen's model and used a three-dimensional (3D) steady Reynolds-averaged Navier-Stokes (RANS) method to analyze the impact of geometrical simplifications on resulting wind loads. The four models for the container ship were selected and step-by-step simplified in geometry from a simple block to the detailed model. Andersen's wind tunnel was numerically modeled by Fluent™, and the effect of blockage in the wind tunnel was investigated. Majidian and Azarsina (2018) also employed CFD for the verification of wind loads according to Andersen's test data. RANS solver ANSYS CFX™ was used, and appropriate conformity of the results was achieved in 10 different simulated loading configurations in head wind. Furthermore, they separately investigated the effect of aft superstructure configuration on air resistance and found that the height of the first bay after the superstructure has the most significant effect.
The reports by Steen (2017a, 2017b) are examples that show the importance of the effect of oblique wind loads on the performance of ship simulators. They introduced metamodels for simulating ship maneuvers. The metamodel is extracted from in-service recorded data of the vessel. However, unknown forces of side winds can increase disturbances and thus increase the uncertainty of the model and simulation results. Among several inputs, wind direction and speed were observed as influential factors to the total uncertainty of maneuvering results.
Previous CFD works on container ship aerodynamic simulations were impaired by one or other of the following: (i) Pressure distribution over the 3D object is not discussed. (ii) Details of numerical model setup are not reported; this particularly includes the boundary condition at the air-sea interface. (iii) Practical importance of the dominant range of oblique winds in a real seaway is not discussed.
With the objective of addressing these drawbacks, the composition of the present research is as follows.
• CFD analysis of an airship. A computer model of Akron airship (Freeman 1932a) is constructed and used to verify numerical results versus US National Advisory Committee for Aeronautics (NACA) test data. Numerical setup is verified to be accurate to the level that it can predict pressure distribution over the surface of the 3D object. This advances the state-of-the-art of using CFD for oblique wind simulations. To the knowledge of the authors, NACA test data is a unique source that contains pressure measurements at nonzero wind angles (Section 2). • Statistical analysis of relative wind directions measured on a voyage of a container ship. The voyage relates to a 2748 TEU fully loaded container ship M/V Azargoun from the port of Durban, South Africa, to Bandar Abbas, Iran. This indicates the practical importance of oblique winds (Section 3). • CFD analyses of a container ship model with different modeling approaches for the air-sea interface. Simulations are within the statistically important range of wind angles as were discussed in a previous section. Results are verified against test data. The innovative approach in modeling the air-sea interface also advances our knowledge (Section 4).

Materials and method
Regarding the present literature and reported technical data, two different methods can be employed for simulation of wind force on a ship and specifically a container ship to date. In the first method reported repeatedly in the previous literature by researchers such as Hassan, White, and Ciortan (2012), Wnek et al. (2015), and Majidian and Azarsina (2018), simulations have been structured based on the presence of a ship in a real environmental condition such as considering the presence of sea surface and opening boundary condition in the top of domain. Meanwhile, in the second method, direct simulation of the wind tunnel is considered (Janssen, Blocken, and van Wijhe (2017) and Taherkhani (2015)). However, the second method takes advantage of the detailed monitoring of pressure on the body compared to the first method, particularly while the experimental data is available for the same simulation that promotes the level of accuracy for further studies.
In the aviation history, airships received the most attention and saw tremendous growth in the early decades of the 20th century. However, in the recent years, owing to their specific properties, increasing trends are observed in the research and design of new airships; e.g., see Trancossi et al. (2015), Zheng and Sun (2016), and Jelenciak et al. (2015).
From 1929 to 1932, a series of experiments were performed in the wind tunnel of the NACA on Akron airship to determine the drag force, lift force, and pitching moment on the bare hull and the hull with appendages (Freeman 1932b). In another set of experiments, a 1/40 scale model of the airship was used to study pressure distributions (Freeman 1932a). The Akron airship model had a length of 5.98 m (19.62 ft.) and had a maximum diameter of 1 m (3.32 ft.). The pressure data were recorded for a nominal air speed of 100 mph, equivalent to 44.7 m/s in the 20-foot (6 m) propeller-research wind tunnel (Freeman 1932a). The pressures due to wind flow were measured on the surface of the bare hull at about 400 locations distributed longitudinally over 26 transverse stations (see Figure 1), at eight pitch angles of 0°, 3°, 6°, 9°, 12°, 15°, 18°, and 20°and two air speeds of about 70 and 100 mph (31.3 and 44.7 m/s). Azarsina et al. (2007) employed modern numerical methods in order to reanalyze the data from the "Akron" pressure experiment.

Mathematical formulation
An unsteady Reynolds-averaged Navier-Stokes (URANS) method was used to solve the governing equations in this study. These mass and momentum conservation equations were solved by commercial CFD software ANSYS CFX™. The averaged continuity and momentum equations for incompressible flows are given in tensor notation and Cartesian coordinates by Equations (1) and (2) (ANSYS 2009).
where ρ is density, u i is the average Cartesian components of the velocity vector, ρ 0 u i 0 u j are the Reynolds stresses, and p is the mean pressure. τ ij are the mean viscous stress tensor components, as shown in where μ is the dynamic viscosity.
In this study, CFD has been employed for computation of the mass conservation, momentum, and energy equations around the body surface. By splitting up (discretizing) the nominated area for analysis into smaller parts and applying boundary conditions, a linear system of differential equations is obtained and solved for the velocity and pressure field. To this end, ANSYS CFX 16.0 is used, which is recognized as one of the most powerful fluid flow analysis software. As an exclusive feature of the software, the program accompanies with the Integrated Computer Engineering and Manufacturing (ICEM), providing one monolith environment. In this way, all steps of creating geometry, meshing, and physical definition of model features, as well as pre and post solver are involved within one unit. A finite volume method (FVM) is used by the software.

Computer model setup
The computer model of the airship is similar to a physical model as described in Freeman (1932a) and Azarsina et al. (2007). The computational domain is designed to be large enough in dimension to minimize the effect of lateral boundaries on the airship. The length of the tunnel duct is 36 m, i.e., two airship lengths forward of its nose and three lengths aft of its tail end. The tunnel diameter is 6 m, similar to the NACA facility (See Table 1). The airship with no pitch angle is aligned on the centerline of the tunnel. The wide range of unstructured 3D grid was generated inside the domain and over the airship surface, including tetrahedrons and wedge cells. Details of mesh quality are presented in Table 2 and Figure 2. Regarding the available features, the numerical setup is summarized in Table 3.
The parameters for thickness of boundary-layer approximation were adapted according to Freeman (1932a). Furthermore, in order to check grid independency, the mesh adaption item (program feature with advanced setting) was activated during running process. Cell size from coarse to fine rendered 7 up to 11 millions of cells.
Consequently, accuracy of results were obtained to a thousandth. The basic SST model accounts for the transport of the turbulent shear stress and gives highly accurate predictions of the onset and the amount of flow separation under adverse pressure gradients. Hence, this turbulence model has been used in this task. Amid, the advection scheme has been set to secondorder high-resolution for convective and viscous terms of the governing equations. In this way the oscillation on the convergence of results are mitigated by choosing blend factor regarding to upwind node behavior (ANSYS 2009). The estimation of wind loads has been performed by the solver of ANSYS CFX 16.0 using a cluster equipped with Intel® Core™ CPU i7-quad core CPU@ 2.50 GHz, 8 GB of RAM. The execution of the analysis of each particular flow simulation took about 8 hours.

Correction for blockage
Generally, in the wind tunnel, the model partly blocked the flow, and this might affect the results substantially; hence, a correction must be considered for decreasing the effect of blockage. International Best Practice Guidelines suggests two terms regarding this issue: first, a maximal blockage ratio (BR = A airship /A wind tunnel ) of 5% for wind tunnel research and 3% for CFD research (Tominga 2008;Franke et al., 2011) and, second, a directional blockage ratio (DBR = Length model /Width wind tunnel ) less than 17% (Blocken 2015) must be considered. These guidelines just partially fulfilled the BR and did not fulfil the DBR in any case (Table 4); thus, a correction must be applied to compensate for the blockage effect.
The force and moment results are corrected for blockage by following the approach of the where C c is the corrected coefficient, C s is the simulated coefficient, m is an expansion factor for the wake, which is extracted from the tabulated data in Engineering Sciences Data Unit (1980), A is the cross section of the wind tunnel (28.26 m 2 ), and S is the projected area perpendicular to the direction of flow, which comprises terms of the frontal projected front area as S x = A f cos α j j and the lateral projected side area as S y = A s sin α j j, where α is the wind angle of attack (AOA). For the turning moment, S v = S x + S y is used. Similar correction factors are used for the container ship simulation, which will be discussed in the next part of this paper, following the values for the projected front and side areas A f and A s of the container ship.

Convergence of results
Aiming to attain fine convergence, with respect to the present software and hardware performance, minimum 10 cycles and maximum 50 cycles of iteration were considered for the simulation of airship, which would stop if the difference between two successive RMS values became less than 0.0001 (see the curves in Figure 3). Since the container ship geometry is relatively complicated compared to the airship, the iteration value was increased to 100 for container ship simulation in Section 4.

Pressure distribution on the airship surface
As mentioned earlier, in order to fulfill the guideline requirement concerning the directional blockage factor, the virtual wind tunnel diameter must be less than 17%, which is basically far from the size of the NACA wind tunnel. It means that the tunnel diameter should extend at least up to 36 m, equivalent to 16.6%, to fulfill the minimum requirement for suppressing the correction factor. However, such a wind tunnel would be quite time-consuming and expensive, particularly for more complicated models. In this respect, one idea is examined to model the airship in a tunnel with 36 m diameter for suppressing any blockage correction automatically and to evaluate the pressure value on the body, drag force, lift, and pitching moment compared to the NACA report directly. The numerical wind tunnel diameter must be reduced to the NACA wind tunnel, that is 6 m, and the results must be compared after applying blockage correction. The obtained results in line with NACA data can lead to a comprehensive insight into the simulation of a model in a virtual wind tunnel practically and provide an optimum solving condition for a container ship in the next part of this paper. The airship is a streamlined body, and separation of flow moves forward as the pitch angle increases. Thus, at inclined angles, the numerical model should be capable of simulating flow separation and large turbulence. To make comparisons, station #9 in the front part of the airship, station #15 in the mid-body, and station #25 in the aft part of the airship are chosen (see Figure 1).
Pressures were measured at different azimuth angles around each station and several pitch angles that are rotated around the center of buoyancy of the airship located 2.755 m from the airship's nose toward the astern (Freeman (1932a), Azarsina et al. (2007)). The azimuth angle varies from 0°at the airship keel to 180°at the top-line.
Five pitch angles were used for simulation: 0°, 3°, 6°, 9°, and 12°. Comparison between numerical results and test data is shown in Figure 4(a-c), which, respectively, correspond to stations 9, 15, and 25 along the airship hull. In the figure legends, e.g., "st9" means station No. 9, the number in subscript is the respective AOA, and extensions "exp" and "num" indicate experimental and numerical results, respectively. It is clearly seen that the trends are similar; at the tail end however, the distance between the numerical and test data becomes larger. This might be because flow turbulence has become too stochastic there and errors are more stressed. The same patterns were observed at other stations along the airship hull. (Relative errors between pressure test data and numerical results are presented in Table A1 of annex.) Furthermore, the drag and lift forces and pitching moments were investigated at eight pitch angles of 0°, 3°, 6°, 9°, 12°, 15°, 18°, and 20°, which show fair conformity with the experimental results, as illustrated in Figure 5(a-c). Drag force from numerical simulation is larger since this is a viscous flow solution including shear stresses, while the test data includes the normal stresses (pressures) integrated over the airship's hull surface area. This difference is larger at smaller pitch angles, where shear effect is more pronounced compared to pressure drag.
Verification of pressure values in the present research promotes the level of accuracy beyond the usual check of only drag and lift coefficients. Next, this numerical setup can be employed for the simulation of a container ship in oblique wind within a small range of AOAs. Before that, however, the practical importance of small ranges of oblique winds has to be studied.

Analyzing a container ship's sea voyage
To use a real scenario in the computational condition, a set of information of one relatively distant sea voyage for a 2748 TEU fully loaded container ship from the port of Durban in South Africa to Bandar Abbas in Iran has been collected and analyzed. The principal information of M/VAzargoun during the voyage is available in Table A2 of the annex.
The literature on this subject is relatively scarce. Among the few, Chen, Shigeaki, and Sasa (2015) showed the significance and safety aspects of constructing a numerical navigation system that can take into effect real environmental conditions to predict a safe route for the vessel. Their numerical navigation scheme was a pilot for a container ship model. Steen (2017a, 2017b) show the justification of a simulation model that is based on full-scale maneuvering trials under the action of wind, wave, current, and steering parameters.
The voyage commenced in Durban at 01:30 LT on 31.05.2017, then completed in Bandar Abbas at 22:30 LT on 10.06.2017, which totally took 11 days. The wind information is extracted from the ship's deck logbook. Normally, this information is mandatory to be written in the deck logbook, which is arranged in the form of a four-hour interval report, unless the ship's duty officer faces any unexpected situation during voyage. The true wind direction must be normally recorded in the logbook in the form of compass arising out of 16 points. Thus, the wind direction is interpreted for 16 compass directions, although wind direction may be slightly different in the actual condition. The wind speed must be logged in the logbook in the form of sea Beaufort scale from 0 to 12 (Moonesun 2010). The dedicated wind speed for each force scale is normally variable from the lower to the upper limit, but for the purpose of this paper, the average speed has been adapted for every force scale, which is presented in Table A3 of the annex.
The ship's course to steer and the respective speed are assumed as the ship's velocity vector, and the wind direction together with the wind force scale is assumed as the wind velocity vector. The resultant vector can be derived from the wind and ship velocity vectors (i.e., relative wind) (see Figure 6).
The resultant wind vector direction has occurred within a range of 0°-40°for 203 h out of 262 total hours of the ship's running time period of the voyage made, which was about 78% of the total running time, and the maximum resultant speed experienced was 24.8 m/s at 39°(See Figure 7). As far as could be observed, the maximum range of wind AOA comprises the range of 0°-40°. Moreover, the sample of Comprehensive Ocean-Atmosphere Data Set (COADS) data (November 1986) showed that over 40% of  merchant ships' observation were obtained for relative wind direction within ±25°of the bow (Moat et al. 2005). Henceforth, using the verified computational domain, the resultant wind vector within 0°-40°r ange was focused in simulation to demonstrate a recurring range of wind AOA on the container ship.
The resultant vector magnitude and direction for all running times represent the dominant direction and magnitude of wind velocity vector on the container ship (see Figures 8 and 9).
These statistics genuinely represent the significance of creating an accurate numerical model for predicting the wind force in quarter head winds, which play an important role in practical wind resistance consideration.

Simulation of a container ship in oblique winds
A fully-loaded post-Panamax container ship with 9000 TEU, which is the same model that was employed in the experimental investigation by Andersen (2012aAndersen ( , 2012b and formerly simulated by Majidian and Azarsina (2018), is used in the present study. The principal dimensions of the container ship and the stacking configuration are presented in Tables 5 and 6. There are several reasons for choosing this post-Panamax container ship, such as access to the precedent experimental and numerical results, wide windage area, and enormous amount of fuel consumption, which is important matter of concern for all involved parties. The container ship is scaled 1/56 to adapt with the domain size that was introduced in Section 2. Several previous works utilized different approaches to model the air-sea interface. All of the earlier studied information regarding the type of domain and the method of air-sea interface are distinctively summarized and presented in Table 7. In this paper, four different approaches are used to model the wind load on a container ship. Figure 10 shows the computer model in four scenarios. Two sets of simulation have been performed based on the wind tunnel experimental test, and two sets of simulation have been performed based on the theoretical approach, with a total of four scenarios, while the airship's solving condition in Section 2 is considered in association with the details of the investigated method.

General conditions of simulation
Grid size is refined especially in the ship's hull, container shape, and geometry because of the presence of sharp corners and edges. In this way, the total number of meshes and nodes increasingly changed and it reached between eight and nine millions for different modules, which were inflated in sharp corners and container surface, whereas the grid independence was carried out up to 11 million meshes for one of the scenarios. The boundary condition also followed the same pattern as for the airship, except for the sea surface, which has been elaborated and explained at each relevant module. The applied meshing and boundary conditions are almost the same for all models except where clearly stated in each part. International Best Practice Guidelines is used to correct the results for effect of flow blockage in the tunnel by using the same procedure as described in Section 2.3 of this paper (Tominga 2008).

Loads and coefficient definition
Wind forces on the post-Panamax container ship including freeboard above waterline, superstructure, aft and fore deck, and loaded containers on deck were simulated; 40 ft container size was used (Shipbusiness 2017). A right-handed Cartesian coordinate system that is fixed across the intersection of the baseline passing through the deck and passing athwartship   0   10   20   30   0130  0900  1700  0100  0900  1700  0120  0900  1700  0100  0900  1700  0100  0900  1700  0100  0900  1700  0100  0900  1700  0100  0900  1700  0100  0900  1700  0100  0900  1700  0110  0915     through the deck at the center of rotation, which is placed halfway the length between perpendiculars, is implemented. To study the effect of wind AOA, the ship is turned counterclockwise in angles of 4°, 8°, 12°, 16°, 20°, 24°, 28°, 32°, 36°, and 40°. The axis directions have been adopted from ITTC 1993 guidelines (Andersen 2013). Moreover, the angle between the ship's advance and head wind is 180°. The axes are defined as follows (see Figure 11): • the x-axis is positive forward, • the y-axis is positive to the port, and • the z-axis is positive upward.
The forces studied in this paper are as follows (see Figure 11): • longitudinal (surge) force X, • lateral (sway) force Y, and • moment about z-axis, yaw moment N.
Drag force is along the relative wind velocity vector and lift force is perpendicular to this direction. The measured forces and moments are postprocessed into nondimensional coefficients to make the results independent of wind velocity and ship size. The following equations have been used with the flow  Figure 11. Indication of the coordinate system. velocity wind profile (U), the density of the air (ρ), ship frontal area (A F ), ship lateral area (A S ), and ship volume (V S ) in a suitable exponent.
Direct representation of the sea surface Numerical setup for airship simulation in a wind tunnel implies free flow condition along the whole wind tunnel; thus, completely the same condition could be considered for container ship simulation. Besides, the study by Majidian and Azarsina (Majidian and Azarsina 2018) represented good conformity with the experimental results for the same container ship with different stacking conditions in the sea-surface interface modeling at 0°AOA. No-slip wall condition has been imposed for sea surface with the velocity equal to the profile's velocity along the axis of the tunnel (see Figure 10(a)). The formation of wind flow around the model is shown in the form of flow streamlines and pressure contours in Figure 12.

Simulation of the model with the image method
According to White (Blendermann 1996), in bounded flow, the effect of distortion is inevitable between the bottom of the object and the lower surface of the boundary. A potential flow solution for this wall proximity is to omit the wall and insert the image of the object as if it is mirrored on the wall. In this respect, as a new innovative numerical trial in a ship's aerodynamic simulation, such a 3D mirrored object was created (see Figure 10(b)), and then results are attained and presented for comparison purpose with the other condition of simulation (see Figure 13). Although this theory is originally applied to potential flows, success has been achieved here with viscous flow conditions in eliminating the undesirable lift force (along the vertical direction); due to the presence of the sea surface, such a lift does not happen in real conditions but a 3D model of the container vessel as is put in the middle of a wind tunneleither experimental or numericalendures such an unrealistic lift force.

Simulation of the model with a rounded base plate
A series of experimental tests were performed by Kim et al. (2015), a group of researchers in Korea Research Institute of Ships and Ocean Engineering, on the container ship with 5000 TEU capacity in the wind tunnel in order to investigate the effect of wind load and possible drag reduction. Subsequently, they compared the test results with CFD (Kim et al. 2015). Regarding this matter, the authors used the same wind tunnel simulation jointly with the NACA wind tunnel setup in order to optimize the domain for the evaluation of axial and lateral forces on the container ship. Similar setup as the Korean base plate was used in wind tunnel simulation. The only difference is that the Korean reference used a rectangular tunnel, while the verified tunnel, with reference to NACA tests, is a round tunnel (see Figure 10(c)). It ought to be noted that the ground plate has been considered in the free slip condition. The indication of domain, velocity streamlines, and pressure contours is shown in Figure 14.

Simulation of the model with a thin turntable
Regarding the real condition of the wind tunnel test, which was performed by Andersen for a 9000 TEU, 1/ 450 scale model (Andersen 2012a(Andersen , 2012b, it was tried to simulate the same wind tunnel setup intertwined with the NACA wind tunnel domain in order to optimize the domain for evaluation of axial and lateral forces on the container ship (see Figure 10(d)). The container ship was placed in the middle of a turntable with the center point located 0.79 m away from the inlet of the measurement section in Andersen's study; however, in this simulation, the location of the turntable is considered two ship's length downstream from the inlet to prevent wind profile distortion. Figure 15 shows the results.

Summary of results
In this section, CFD simulation results of the four scenarios -flat surface, image method, rounded base plate, and turntablefor drag force, side force, and yawing moment coefficients are compared with experimental data from Williams et al. (2006), Blendermann (Maciejewski and Osmólski 2002) and Kim et al. (2015). As can be observed, simulation results for scenarios 1 and 3 are very close to each other. The reason is that both scenarios model the sea surface as a stationary extended plane; the only difference is that scenario 3, which mimicked the wind tunnel test in Kim et al. (2015), had a rounded front. Scenario 2 which is the innovative image method shows quiet close similarity with the Andersen's test results. Results from scenario 4 show variations that might stem from the effect of attached turntable base beneath the model.
In comparison to the test data from four different sources, it is observed in Figure 16(a-c), that scenarios 1 and 3, where the sea surface was represented by a flat extended plane, overestimate the experimental results, while scenario 2 shows higher conformity to the experimental results. One may question the difference in results from scenario 3 and its experimental counterpart (Kim et al. 2015). The reason may stem from surface roughness and turbulence phenomena, which are not identical between the physical and numerical models. On the other hand, simulation of a ship on a turntable resembles the mounted outfit in the actual experimental test in a wind tunnel, seemed not to be a justified hypothesis as the results stand quiet far from other scenarios/experimental data especially in axial force where its curve leaped drastically and in the yaw moment when it partially converged.
Scenario 2, an innovative approach to model a surface vessel inside a 3D wind tunnel, is selected as the most effective model. The results are in close agreement with the experimental data reported by Andersen (2012b) as indicated in percentage errors in Table A4 in the annex. For the yawing moment in Figure 16(c), the test results of Kim et al. (2015) are not in agreement with those from Figure 12. Indication of velocity streamlines and pressure contours around the model at 12°"AOA." Andersen (2012b), and CFD results are also comparable to those of Andersen (2012b). The reason might have been that containers' spacing and details of the experimental setup were in a way that although the axial and lateral forces are similar, the yawing moments are different. In other words, center of effortaxial location of an axis about which the pitching moment is zerois different in two physical models (Williams et al. 2006). The authors did not have access to such details in performing experiments, and Andersen (2012b) has been taken as the reference to evaluate simulation errors in Annex A.
To a certain extent, the double-model method can take lead among all of the tested scenarios here and the previous literature for evaluation of wind load on ships in a virtual wind tunnel. The reason is that the double-model is rarely affected by the presence of appendage in the tunnel and prevents the dissipation of wind profile at height or over the flat plates that are commonly embedded as a sea surface. Furthermore, compared to the presentation of a single model in the middle of a tunnel, the undesirable lift force is omitted. The downside of this scenario is that twice as many mesh elements are produced over the vessel surface, but many more elements are reduced as a flat plate (scenarios 1 and 3) or a base (scenario 4) is removed. Here, a drag coefficient is also meaningful to be evaluated for the aerodynamic performance of the container vessel. A drag coefficient is defined as follows: Results of scenario 2 as reported above were used and substituted into Equation (8). The resulting drag coefficient was fitted with a sum of two sine functions as follows (see Figure 17): C D ¼ 15:84 Â sinð0:04006 Â α þ 1:016Þ þ 15:06 Â sinð0:04207 Â α þ 4:189Þ where wind angle α is substituted in degrees, and the resulting terms in parentheses are in radians. Use of sum of sine functions has the advantage that the fitted curve changes a curvature: by increasing the wind angle, the drag coefficient cannot increase unbounded. Equation (9) was derived in a range of 0°-40°, where goodness of fit is measured as SSE: 0.004491, R-square: 0.9983, adjusted R-square: 0.9967, and RMSE: 0.02997. Respective data for the drag coefficient have been listed in Table 8. Results for this range of AOAs have already been presented in Section 4.7, Figure. 16(a-c). In the section Analyzing a Container Ship's Voyage Data, the surge force indicates the additional aerodynamic resistance due to wind flow. Here, the authors suffice to summarize an added aerodynamic resistance for where A XV is the area of maximum transverse section exposed to the wind, C X is the wind resistance coefficient, V WR is the relative wind speed, ρ A is the mass density of air, and ψ WR is the relative wind direction. Aerodynamic resistance versus wind angle is shown in Figure 18.

Summary
The significance of this research can be observed from the following factors.
• A numerical wind tunnel has been verified to predict pressure distribution over a 3D object. The NACA test data for Akron airship could be predicted with acceptable accuracy up to 20°pitch angle. • Log data of a container ship's voyage was analyzed and the range of oblique winds from 0°to 40°was concluded to be statistically significant in aerodynamic loads.   Table A4 in the annex).
• Four scenarios to interpret the effect of air-sea interface, inside the previously verified wind tunnel, were attempted. They respectively included: (i) a rigid bottom surface of the tunnel, (ii) a mirrored 3D object of the vessel mounted on the tunnel centerline, (iii) a base plate with a rounded front, and (iv) a circular turntable. Of course, modeling the effect of sea surface is an important task in this problem. Ultimately, the image method, as an innovative method of the present study, provides more accurate results compared to the other simulation schemes, and this can be introduced as a new method for the same task in further CFD studies.   Figure 18. Aerodynamics-added resistance versus wind angle.
• Numerical simulations for extreme cases of the voyage data were performed.
In future, the authors would conduct a physical model test. In addition to force measurements, pressure transducers can record the pressure distribution over the scale model of the container ship. ANSYS CFX provides graphical outputs for vortices and the wake flow. The authors have examined this feature and are working on the quality of wind field around a container ship while it maneuvers in a surgesway-yaw plane.
Cargo configuration effect on oblique wind loads is another important by-product of the present research.
Last, but not least, is an interactive model of waterfree surface, which can respond to wind shear stresses. Such a model requires supercomputing capacity and consequently will add to the accuracy of the numerical model. Table A1. Relative errors (percentage) between pressure test data and numerical results at several azimuth angles around stations 9, 15, and 25 (graphical data are shown in Figure 4).