Mesh Properties for RANS Simulations of Airfoil-Shaped Proﬁles: A Case Study of Rudder Hydrodynamics

: A good mesh is a prerequisite for achieving reliable results from Computational Fluid Dynamics (CFD) calculations. Mesh properties include mesh types, computational domain sizes, and node distributions. However, in literature, we found no clear consensus about what these properties should be. In this article, we performed a case study on ship rudders to determine what the suitable mesh properties are for airfoil-shaped proﬁles. A classic NACA 0012 proﬁle is chosen as an example, and commercial packages ANSYS ICEM are applied for meshing with an ANSYS Fluent solver. With a strategy in consideration of relationships among different mesh properties, a comprehensive parametric investigation is conducted to study the impacts of these properties on the accuracy of rudder hydrodynamic coefﬁcients obtained by CFD methods. The step-by-step study outputs recommended Reynolds numbers, domain sizes, and near- and far-ﬁeld node distributions for mesh types with distinct topology structures, i.e., C-mesh, O-mesh, H-mesh, and Hybrid-mesh. Speciﬁcally, the study shows that a critical Reynolds number is needed for the perspective of efﬁciency, while a domain extending 60 times of the chord length enables the boundary effects to be negligible. As for node distributions, the near-ﬁeld nodes should be treated carefully, compared with those in the far-ﬁeld. After that, corresponding mesh properties for different calculation objectives are illustrated in detail based on the characteristics of mesh types mentioned above. With the proposed strategy for mesh reﬁnements, impacts of different mesh properties on rudder hydrodynamics are clariﬁed and recommended settings are applicable for other airfoil-shaped proﬁles such as wind turbines and marine propellers.


Introduction
Stimulated by developments in computer power and mechanics theories, Computational Fluid Dynamic (CFD) methods were widely applied in ocean engineering to study the hydrodynamic performance of marine structures. Generally, computation domains containing target structures need to be discretized into multiple mesh elements to solve differential equations numerically. In ship maneuvering simulations, the presence of CFDbased tools enables researchers to obtain more precise hydrodynamic coefficients compared with that of traditional empirical methods, which is beneficial to establish individualized and accurate mathematical maneuverability models for different objectives. As the most widely used steering devices for ships, rudders play an important role in the maneuvering of ships [1], and the accuracy of lift coefficients has a great impact on motion predictions. Many CFD simulations were conducted to study rudder hydrodynamics (Liu et al. [2], Badoe et al. [3], Van Nguyen and Ikeda [4]), but few concentrate on mesh generations for rudders. In this paper, the study on mesh properties for RANS simulations of rudder hydrodynamics is performed to present detailed insight into meshes around rudders.
A good mesh is fundamental to the convergence and the accuracy of solutions, while a poor one leads to nonconverging results, inaccurate solutions, and an increase in computation time. Due to the complexities of geometries in structural surfaces and the lack of automatic generation methods, it costs plenty of time to generate meshes in calculation processes. On the one hand, mesh quality is influenced by various factors, but on the other hand, it is unrealistic to adjust all of them to the best states in practice, which requires more computation power and computation time.
A good quality mesh should have a sufficiently large number of cells, implying that a further increase in the number of cells will not lead to significant changes in the solutions. The mesh needs to be fine enough in the regions where fluid characteristics change rapidly. Therefore, mesh properties need to be coordinated to make a trade-off between the accuracy and the time scale. In literature, different mesh properties, i.e., mesh types, domain sizes, and node distributions, are applied in RANS simulations of airfoil-shaped profiles. Lutton [5] indicated that the O-mesh is superior to the C-mesh in determining the pressure coefficients in the vicinity of leading and trailing edges, while the C-mesh shows a better resolution in the wake. Basha [6] constructed a structured O-mesh, a structured C-mesh, and a hybrid mesh on a NACA 0012 airfoil, while the hybrid mesh shows better results in drag coefficient predictions due to finer mesh resolutions in the boundary layer domain. Stuck et al. [7] selected a hybrid mesh with quadrilateral unstructured cells rather than triangular ones to achieve smoother transitions from the boundary mesh to the unstructured outer mesh for RANS simulations on a rudder profile. Parametric grid independence studies are conducted by varying chord-wise and layer-wise grids, respectively, and optimum values are chosen based on variations of integral force coefficients. Turnock et al. [8] carried out decoupled studies on the boundary location and mesh node distributions in different directions for a NACA 0012 section, while related parameters are determined by convergence study.
As no clear consensus about mesh properties was found in current studies, systematic parametric investigations on mesh properties are performed in this paper. As naval architects, our primary interest in CFD lies in applications related to the flow around ships and resultant hydrodynamic forces. In this particular case, we focus on the flow around the rudder and the resulting rudder performance. Kim et al. [9] stated that the NACA series are the most widely used economic rudder profiles. Thus, a classic NACA 0012 airfoil-shape rudder profile is utilized to study mesh properties for its simplicity in the shape, which makes mesh properties easier to change and enhances the applicability of obtained results for other complicated structures.
This paper analyzes the impacts of mesh properties on RANS simulations and proposes suitable meshing strategies for rudder hydrodynamics. Following this introduction, Section 2 further discusses the reason why a 2D rudder profile is selected as the research object in this paper. Section 3 explains definitions of mesh properties in this study. Section 4 illustrates applied research methods in detail, and Section 5 discusses the impacts of different mesh properties based on simulated cases. Finally, a recommended meshing strategy is presented to provide guidance for hydrodynamic analysis in Section 7.

The Selection of the Research Object
Bare hulls, propellers, rudders, or full-appendages ships were all studied in CFD methods in literature, though not all of them are suitable for mesh properties studies. There are great curvature changes on the surface of hulls and propellers, and that is why mixture mesh topologies and millions of cells are desired to capture their flow features. The complex nature makes it difficult to perform parametric investigations on a single mesh property economically. What's more, the geometries of hulls or propellers in the marine industry vary from each other significantly, and rules derived from limited cases may not be applicable for those with different configurations.
Although special rudders drew a growing interest, conventional spade rudders, or those modified but still geometrically similar versions like flap rudders, are still mainstream maneuvering devices, however. Besides, the meshing strategy aiming at a NACA profile can be easily applied to its geometrically similar counterparts. Force coefficients of rudders during ship maneuverings can be obtained based on 2D profiles in open water conditions, which also makes the investigation of 2D profiles beneficial to maneuverability studies. For all these reasons, in this paper, a 2D profile, in other words, a rudder with an infinite aspect ratio, is determined as the research object. According to Fujii [10], Fujii and Tsuda [11,12], the normal force coefficient of actual rudders whose aspect ratios are limited can be expressed as: (1) where C N is the rudder normal force coefficient, and 6.13 is an approximate constant to estimate the partial derivative of C N to sin α R of 2D sections. Liu et al. [2] calculated ∂C N ∂ sin α R for different rudder profiles by performing 2D CFD simulations, which can modify the original transformation relation in Equation (1) and obtain force coefficients for 3D rudders in open water conditions. Since rudders generally operate in the propeller slipstream, rudder hydrodynamics affected by the propeller can be estimated based on series model tests by Nienhuis [13], and the experimental data are shown in Figure 1. Compared with that of open water conditions, the stall angle for the rudder increased to about 35 • with propeller impacts, while ∂C N ∂ sin α R was not significantly changed. Hence, based on 2D simulation results, corrected formulations for a rudder with an aspect ratio Λ G considering propeller impacts can be expressed as: D 0 can be obtained from 2D CFD simulations, and k p can be estimated experimentally. The availability of Equation (2) in ship maneuvering studies was validated by Liu et al. [2,14] against two typical inland vessels. Apart from the transformation relation, another factor that makes the 2D study effective is that spanwise mesh generations can be omitted and cells in the profile section can be focused on. Therefore, 2D rudder profiles are set as the research object in following studies.

Mesh Properties
In literature, mesh types and computation domain sizes are generally determined before studying grid independence on certain mesh parameters, and some settings are shown in Tables 1 and 2. Mesh properties in grid independence study are denoted as near-field and far-field parameters according to locations of distributions. All of these mesh properties are studied in this paper. The detailed descriptions of mesh properties above are in Sections 3.1-3.3.

Mesh Types
Based on the connectivity of elements, meshes are identified as structured, unstructured, and hybrid. As regular connectivity can be expressed as a two/three-dimensional array, the structured model is highly efficient in the usage of computational resources. However, regular connectivity restricts the element type to quadrilateral in 2D and hexahedral in 3D, which inherently limits applications of the structured mesh for complex geometries. For simple geometries, a structured mesh may have better convergence and higher resolution than an unstructured mesh.
Instead of a single-block structured mesh, a block-structured or multiple-structured mesh is more widely applied. Three commonly used topologies, C-mesh, H-mesh, and Omesh, are shown in Figure 2. It is superior in computational efficiency to the unstructured mesh and more flexible in handling the complex geometries than the single-block structured mesh. The choice of the mesh topology has a considerable influence on mesh quality [22], which depends on the domain geometry, the structure of the solution, and the topology in the adjacent block. According to Liseikin [22], in H-mesh, the computational domain is square and two singularities exist in the interior boundary. O-mesh and C-mesh are both in solid squares, while the interior boundary for the former one is smooth and the latter has a singularity the same type as the H-mesh. In this particular case, some modifications are made based on basic topologies above to accurately capture the characteristics of the flow near the rudder, as shown in Figure 2. For brevity, modified types are still marked as C-mesh, H-mesh, and O-mesh, respectively. The irregular connection of unstructured meshes allows for more freedom in element choices, typically triangular in 2D and tetrahedral in 3D. Compared to that of the structured mesh, the unstructured mesh is more suitable for complex configurations, but it can be highly inefficient. Aftosmis et al. [23] reported that unstructured triangular meshes are 50 times more expensive in both memory and time than that of structured quadrilateral meshes with nearly the same accuracy. Unstructured meshes, however, have advantages over structured meshes in mesh refinement and generation time. In practice, hybrid meshes are applied instead of pure unstructured meshes. Hybrid meshes commonly consist of a structured portion near wall surfaces to capture the boundary layer and an unstructured portion that fills the domain. Thus, hybrid meshes inherit the advantages of both structured and unstructured meshes, such as good orthogonality to wall surfaces, suitability for mesh refinement, flexibility for complex geometries, and fast generation. Figure 3 demonstrates a hybrid triangle mesh, which has quadrilateral elements for the boundary layer and triangular elements for the remainder of the domain.

Domain Sizes
The computational domain is the geometrical region that bounds the numerical simulation [17]. To obtain highly accurate solutions, the position of boundaries has to be discussed to demonstrate that the interior flow field is not disturbed. Thus, there is the need to study the domain size, which not only minimizes the influences of the boundaries, but also prevents an excessively large domain. Typical mesh types correspond to related topology structures, which lead to different shapes of calculation domains. The geometric parameters for different mesh types are shown in Figure 4, while subscripts, i.e., in, out, and side, stand for distances from the rudder to the inlet, outlet, and side boundary of corresponding mesh types, and O r is the radius of the domain of O-mesh.

Node Distributions
For structured meshes, since cell distributions vary with node distributions in different directions of the computation domain, convergence studies need to be done by varying related mesh parameters. In this study, node distributions are divided into two categories, including near-and far-field distributions.
Near-field distributions are located in the boundary layer regions near the wall, which are generally refined to capture the flow around the structure. The chord-wise spacing is the distance between two nodes along with the profile, and the layer-wise growth rate is the distances between adjacent mesh points along a mesh line, which determine the boundary layer mesh along with the first mesh height.
Different mesh types share similar near-wall distributions, while it is difficult to define far-field distributions owing to differences in mesh topologies. In this paper, far-field distributions are defined as node distributions that do not contain near-field distributions.

Research Approach
A series of parametric studies are conducted to observe the impacts of mesh properties on the prediction of the accuracy of solutions. The classic validation profile NACA 0012 is tested through commercial CFD code ANSYS Fluent. Strategies for investigating mesh properties mentioned in Section 4.1. Numerical methods utilized to study the hydrodynamic performance of the rudder profile are described in Section 4.2, and the error calculation method is mentioned in Section 4.3.

Strategies
Different mesh properties contribute to mesh quality in varying degrees, and it is the combination of these properties that determines the accuracy of calculation results. Instead of rough global refinement of cells [24], several main impact factors on mesh fineness and the required number of cells are illustrated with simulations. This parametric mesh independence study is chosen because optimization on individual parameters specifically addresses the crucial impact factor of mesh quality, avoiding the waste of cells [7]. The global refinement is more frequently applied for unstructured meshes than structured meshes. However, due to the lack of specific classification on mesh properties, mesh properties above are not independent of each other strictly. Cases with different mesh types also differ from each other in domain sizes, node distributions, and numbers of mesh elements, etc.
As shown in Figure 5, the typical mesh type corresponds to a related topology structure, which leads to different shapes of calculation domains and connection patterns between the profile and domain boundaries. With domain size varying, node distribution in different directions changes. Further, the number of mesh elements is concerned with all mesh properties above.  The strategy in Figure 6 is applied to investigate how different mesh properties affect the calculation results of rudder hydrodynamics. Firstly, cases of C-mesh type with Reynolds numbers in a wide range are tested to obtain the critical Reynolds number (Re critical ), above which hydrodynamic coefficients tend to be stable. Secondly, the domain sizes of the C-mesh are varied to obtain the proper boundary sizes that are large enough to eliminate boundary effects for the following simulations. Thirdly, cases are divided into 4 groups classified by 4 different mesh types. Since it is inappropriate to compare mesh qualities of different mesh types when other mesh properties are not determined yet, the impacts of node distribution on each mesh type should be studied independently. In this part, the node distribution is the only factor that affects mesh quality for the typical mesh type, and the effects of node distributions are studied based on Re critical and the domain size obtained in the previous stage. Therefore, grid independence studies are conducted on parameters in the near-field, i.e., chord-wise spacings and layer-wise growth rates, and the far-field, i.e., radial nodes, wake nodes, and element sizes. Finally, four mesh types with optimal mesh properties are compared to give recommended mesh settings for rudder hydrodynamics.

Numerical Methods
The RANS methods presented in this paper use a time-average Reynolds decomposition, which assumes that all the components of the flow velocity and pressure consist of a mean value and a bounded fluctuation to represent turbulence. The governing equations for incompressible viscous flow are as follows: whereū i is the mean velocity component in the x direction, ρ and µ are the density and viscosity coefficient of the fluid,P is the mean pressure, f i is the body force, while ρu i u j is the Reynolds stress. In RANS simulations, turbulence models are introduced to close equations above. The k-ω SST model was designed to give highly accurate predictions of the onset and the amount of flow separation, which provides a better description of flow around the rudder. Therefore, the k-ω SST model is applied for simulations in the following parts. Numerical settings for ANSYS Fluent is shown in Table 3.

Error Calculation
Routinely, nondimensional coefficients are used to compare the rudder hydrodynamic performance with that of various design choices, while main hydrodynamic characteristics are the lift coefficient C L , the drag coefficient C D , and the moment coefficient C M ([1]), which are expressed as: where L R , D R and M R are the lift force, drag force, and moment, respectively, ρ is water density, α is the angle of attack, V R is the rudder inflow speed, A R is the rudder area, and c is the chord length of the rudder profile. The center of pressure is set as 0.25 c from the rudder's leading edge, based on engineering practice and experimental data [25].
To evaluate the mesh quality, wind tunnel results of the NACA 0012 profile, which are tested by Ladson [25], are taken as the benchmark case. Moreover, three independent compressible CFD codes [19], i.e., CFL3D (NASA LaRC, USA), FUN3D (NASA LaRC, USA), and NTS (NTS, Russia), are cited to make a crosswise comparison with the Fluent solutions. Results of the benchmark and the CFD codes are achieved in essentially incompressible air with a Mach number of 0.15. Thus, they are comparable with the results obtained in incompressible water in this paper. Presented results are compared to that of the benchmark cases in relative differences under various angles of attack at a critical Reynolds number as follows: where ∆C α is the relative differences of the calculated lift, drag, or moment coefficients achieved by Fluent and provided by the benchmark data, respectively, and ∆C(AVE) is the average relative differences for a series of α.

Uncertainty Estimation
The uncertainty estimation is important to evaluate the reliability of data in CFD calculations. According to International Towing Tank Conference [26], simulation numerical uncertainty U SN and simulation modeling uncertainty U SM need to be assessed respectively. In this paper, only the numerical errors are considered since modeling errors are difficult to quantify in practical applications [27]. Generally, numerical errors relate to grid size, time step, and other parameters. This paper performs steady simulations and focuses on mesh properties. Accordingly, only the grid uncertainty U G is calculated for uncertainty estimations. With validation uncertainty U V , which is simplified as U G in this study, obtained, U V and the comparison error |E| can be compared. If |E| < U V , the validation is achieved at U V level.
The Grid Convergence Index (GCI) method [28,29] and the correction factor (CF) based method [30] remain popular in uncertainty estimations, however, it is deficient that both methods fail to give reliable estimations when the order of grid convergence p is too high (i.e., p > 2) or too small (i.e., p < 0.5). In other words, GCI and CF methods require data in the so-called asymptotic range, which means that data on grids are fine enough to give a single dominant term in a power series expansion of the error [31]. Considering that practical flow phenomenons are generally complex, and it is inevitable that p may be not in the range of [0.5, 2], Eça and Hoekstra [31] proposed another GCI method based on least-square roots (LSR-GCI). Such a method allows for error estimations even if the monotonic convergence is not achieved, and it is proven to be more reliable compared with GCI and CF based methods by De Luca et al. [32]. Therefore, the LSR-GCI method [31] is applied for uncertainty estimations in this paper.
For four solutions φ i , (i = 1, 2, 3, 4) obtained from systematically refined cases and characterized by typical cell sizes h i , (i = 1, 2, 3, 4), the discretization error can be estimated with one of the four following equations: where φ 0 is the estimate of the exact solution, α is an undetermined constant, and p is the observed order of grid convergence. Besides, the following functions are proposed to select the error estimator: Details in selecting the error estimator are shown in Figure 7 according to Eça and Hoekstra [31]. φ * 0,RE , α * RE , p * RE are variables that make S RE (φ 0 , α, p) minimum, which can be denoted as Corresponding standard deviations can be expressed as: For no-weighted cases, nw i = 1. Following criteria are applied in selecting the error estimators: Figure 7) Select p * corresponding to min σ * RE , σ w *

RE
• Criterion 2 (S 2 , S w 2 in Figure 7) When p * > 2, small uncertainty quantities can be obtained due to overestimates of the order of accuracy. As recommended by Eça and Hoekstra [31], when p * > 2, the error estimator is chosen by comparing the standard deviations resulting from Equations (16)- (19). Then, ε φ is calculated with the selected estimator. This procedure substitutes the estimated order of accuracy with 1 or 2 because the exponential power of error estimators to be selected (Equations (11) and (12)) is either 1 or 2. Slightly different from Eça and Hoekstra [31], De Luca et al. [32] substitute p as the theoretical order of accuracy p th , i.e., p = p th = 2, when p * > 2. In fact, the method by De Luca et al. [32] is more conservative, and it is applied in following calculations. After that, the error estimator can be selected according to min Figure 7) Select the error estimator according to min σ * , * Criterion 1 * * * > 0, * > 0 , * Criterion 1 * * * < 0, * > 0 * = * * > 0, * < 0 * = * * < 0, * < 0

Results and Discussions
With the strategies and numerical methods illustrated in Section 4, systematic simulations are conducted to study the impacts of mesh properties on hydrodynamic performance of the NACA 0012 rudder profile in this section. For each case with different mesh settings, errors with respect to benchmark data [19,25] are displayed in tabular form. Besides, uncertainty estimations are conducted for cases with varying cell sizes. For colored error tables in this section, the color legends are marked based on the values of relative differences. The smallest difference is the darkest green and the largest one is the darkest red for each mesh group, while other colors between the darkest green to red show the transition. The impacts of Reynolds numbers are presented in Section 5.1, and Section 5.2 discusses the selection of domain sizes taking the C-mesh as an example, and impacts of node distributions and mesh types are discussed in Sections 5.3 and 5.4.

Impacts of Reynolds Numbers
Before conducting parametric investigations, the Reynolds number should be checked first. The Reynolds number (Re) is defined as the ratio of inertial forces to viscous forces, which represents the viscous similarity among flow patterns and can be expressed as: where U ∞ is the inflow velocity, L is the characteristic length, and µ is the viscosity coefficient. For rudder hydrodynamics investigations, due to the limited model size and capacity of the test facility, Re of tests in wind tunnels and towing tanks has to be scaled. Analogously, model-scale rudders are preferred by CFD methods since a fine mesh and a large domain are needed for full-scale simulations. In most cases, the similarity law of viscosity and gravity can not be satisfied simultaneously. Generally, experiments or numerical simulations are conducted under a Re which is large enough to guarantee that the turbulence develops sufficiently. A high Re ensures a fully turbulent flow and reduces the effect of the laminar-turbulent transition, achieving an increase of prediction accuracy compared to that of an analysis at low Re with the same mesh. Cases with larger Re show less calculation efficiency seeing that more cells around boundaries are required. The height of the first row of cells ∆s varies with Re, as shown in Figure 8, while increasing Re leads to decreasing heights. When discretizing domains into multiples mesh elements, the aspect ratio for a single element should be around 1 to ensure mesh quality. When it comes to cells near walls, aspect ratios are connected with ∆s and chord-wise spacings ∆c. As shown in Figure 9, more nodes along the wall are distributed with smaller ∆s to keep aspect ratios compliant with the requirement of solvers. Thus, a decrease in Re leads to a decrease in the number of cells. To sum up, we suggest that using a minimum Re which is achievable in experiments, ensures full-turbulent conditions, and can be solved with fewer mesh elements.

Δs
(m) To determine the threshold of the Reynolds number, above which the turbulence develops completely and force coefficients do not change significantly, cases with varying Re are simulated in this section. Based on the flat-plate boundary layer theory and the calculation method in [33], the height of meshes' first layer ∆s above the rudder surface depends on y + and Re. According to the application scope of the k-ω SST model, y + is set as 1 for all cases, and corresponding heights are shown in Figure 8.
Common benchmark wind tunnel tests are carried out at Re in the range of 1.00 × 10 5 [34] to 1.00 × 10 7 [25]. As low-Reynolds-number RANS analysis is still challenging [35,36] and high-Reynolds-number simulations may be expensive in computation time, the present work tests the NACA 0012 profile at Re in a range of 2.00 × 10 5 ∼1.89 × 10 7 . Among these aerodynamic experimental results, limited cases with Re = 3.94 × 10 6 ∼1.89 × 10 7 and a small Mach number can be taken as validation results of incompressible water simulations, as the compressibility effects of a fluid with a Mach number smaller than 0.20 is small. Table 4 presents the case settings and results of various Re under distinct α series, as specific values of α are listed in [25]. Figure 10 presents calculation results of hydrodynamic coefficients versus the same rudder angles for different Re, which indicates that the lift curve rises with an increase of the Re while the drag curve decreases. Compared to that of the lift coefficients, the drag coefficients are more sensitive to changes in Re. The drag coefficient under 15 • at Re = 1.89 × 10 7 is about half of the value at Re = 4.00 × 10 5 , whereas the lift coefficient is 1.20 times larger. The differences in lift and drag forces between low and high Re increase with an increasing angle of attack. Consistent with findings by Ladson [25] and Molland and Turnock [37], Re = 6.00 × 10 6 can be considered as a threshold value for the mesh generation, above which little variation may be found. Moreover, the y + along the chord for four cases under small and large rudder angles with varying Re are shown in Figure 11. y + values for the case with Re > 5.97 × 10 6 are roughly between 0 and 1, indicating that the k-ω SST turbulence model is applicable for the case. Since ship rudders in propeller slipstream usually experience full turbulence with large Re, low Re conditions are not of interest in current studies. In this paper, we set Re of the following simulations as 6.00 × 10 6 . Table 4. Comparison of results of hydrodynamic coefficients under different Re with that of experimental benchmark data adopted from [25] in relative differences.

Impacts of Domain Sizes
A C-mesh is applied to investigate the effects of the domain size on RANS solutions by conducting simulations with various sizes against varying cell numbers in this section. Test parameters and related relative differences are given in Table 5. C in and C out are doubled for each case compared with that of the previous one, while the narrowest domain extends 7.5 c upstream and 15 c downstream, which is 16 times smaller than the largest ones, i.e., 120 c and 240 c. Since equivalent mesh nodes correspond to different spacing variations in a variety of domains, different nodes distributions are applied for each domain to determine reasonable meshing strategies.
Pressure distributions in the whole domain and around the rudder are presented in Figure 12. From Figure 12a-c, regions with great gradient changes near the rudder can be observed even for the smallest domain. Comparing Figure 12d-f, neither high-pressure regions near the stagnation points nor regions near the trailing edges show significant distinctions. The impact of domain size on pressure distributions around the rudder is relatively small, which corresponds to minor differences of C L mainly derived from pressure integration.
Compared with that of C L , accuracy of C D is more sensitive to domain sizes in Table 5. With an increase of the domain size, the accuracy of the prediction of C D is more significantly improved than C L . Small domain sizes show poor accuracy in C D but obtained C L is acceptable for engineering applications. The relative differences of C D for α = 15.20 • drop from 55.19% to 34.79%, indicating that boundary effects can evidently change viscosity components around the rudder, which mainly contribute to drag forces. On the other hand, larger domains, like 60 c/120 c and 120 c/240 c, tend to show high and stable precision with more cells. However, further expansion of the domain size may dramatically increase the number of cells and computation time. Table 5. Comparison of results of C-mesh based on various domain sizes with that of experimental benchmark data adopted from [25] in relative differences. From the perspective of efficiency, a domain size that eliminates boundary effects and meanwhile requires fewer mesh nodes is desired for following investigations. Table 5 shows that a domain size of C in = 30 c and C out = 60 c with about 2.50 × 10 5 cells achieves a good balance in both accuracy and efficiency, and Table 6 shows that the validation process is achieved for C L obtained from the case (φ 2 ). This domain size is large enough to obtain sufficiently accurate estimates of lift and drag, and the number of cells is still manageable by common desktop computers. Thus, it is applied for the remainder of the simulations in this paper. Table 5 indicates that even more accurate predictions can be obtained with larger domains. Thus, we suggest applying a large domain when computation power and time are available. If the lift is the only purpose, a small domain size of 15 c/30 c may be more favorable as it requires less computation power and time than a large domain.

Impacts of Node Distributions
After determining the test Re and the domain size, raw meshes are ready for the mesh independence study for different mesh types. The common procedure is to carry out simulations on an initial mesh with a residual error in a range of 1.00 × 10 −4 ∼1.00 × 10 −5 . After that, refine the mesh globally to around 2 or √ 2 times the initial mesh. Next, run simulations with the refined mesh and compare solutions obtained from the coarse mesh and the fine mesh. Repeat the refinement until the results do not significantly change with a finer mesh. The mesh independence study requires at least 3 solutions to evaluate the convergence of certain inputs [26]. Considering the computation time, it is always better to use the smallest number of cells. In this paper, mesh resolutions are changed by conducting local rather than global refinements, while the grid refinement factor R G is still around √ 2 for target parameters.
Since the flow near boundaries of the rudder profile changes more violently than that in the far-field and has larger effects on rudder hydrodynamics, near-field node distributions are studied by fixing far-field node distributions for four mesh types. After that, the impacts of far-field distributions are studied with proper near-field distributions.
In this section, with a domain size of 30 c/60 c, impacts of near-field node distributions are investigated for four mesh types, as shown in Figures 13-16. Meshes in different regions are categorized into two parts, which are located in the far-field and near the rudder surface. Near-field meshes are defined as those located along the rudder chord and within the boundary layers regions. The domain for C-mesh in Figure 13 consists of a semicircle and a rectangle, and nodes distributed parallel to the radius are defined as in radial directions, and the rudder in O-mesh is surrounded by radial nodes in Figure 15. The nodes extending from the trailing edge to the outlet of the boundary are defined as in wake directions. Patterned after the definitions above, node distributions for H-mesh are shown in Figure 14, and nodes from the leading edge to the inlet are in inlet directions. For Hybrid-mesh in Figure 16, considering quite different mesh generation methods compared with that of structured meshes, meshes in the refinement area near the rudder are treated as near-field meshes while others are in far-field.

Impacts of Near-Field Node Distributions
In this part, impacts of near-field node distributions are studied by changing chordwise spacings and layer-wise growth rates for C-mesh, H-mesh, O-mesh, and Hybrid-mesh. The near-chord refinement is performed by changing the number of nodes in the chordwise direction for structured meshes and chord-wise spacings for Hybrid-mesh with R G around √ 2. For structured meshes, a reference setting is coarsened or refined to perform parametric studies. The quantitative relation between a coarsened or refined case and the reference case can be expressed as: where i denotes the level of refinement with respect to the reference case, and N c i is the number of nodes along the chord. Taking C-mesh for example, the rudder is wrapped by denser cells in chord-wise directions, as shown in Figure 13. As for cells in the boundary layer region of the rudder, much attention should be paid to capture violently changing velocity gradients near the wall. The layer-wise growth rate is investigated to provide a detailed view of mesh resolutions in boundary layer regions. Cells around the rudder with varying chord-wise spacings for C-mesh are shown in Figure 17. As the chord-wise spacing decreases in Table 7, the number of cells with large aspect ratios in the wake increases, which introduces the difficulty of convergence, especially under large angles of attack. The accuracy of C D is slightly improved with denser chord-wise node distributions, while variations of C L are even smaller. Larger deviations are obtained for C M for small α, probably due to measurement errors in experiments induced by small absolute values. For α = 4.18 • and α = 8.22 • , better mesh resolution along the rudder lead to lower precision, which are related to convergence problems due to the fact that small flow features are being solved. Since calculation results tend to stabilize without diffusing after log R G N c i0 = 2 and 5.48 × 10 5 mesh quadrilaterals show relatively high accuracy, the case with log R G N c i0 = 3 is selected for the following simulations on C-mesh, and corresponded uncertainty estimations are shown in Table 8. Table 7. Comparison of results of C-mesh based on various chord-wise spacings with that of experimental benchmark data adopted from [25] in relative differences with 670 nodes in chord-wise direction for reference setting.     For H-mesh, besides chord-wise nodes, nodes distributed along with radial directions in the refinement region also count in near-field node distributions. Therefore, chord-wise nodes and radial nodes near the rudder are both refined for H-mesh. Due to discrepancies in refinement regions in Table 9, the NOC of H has a more rapid growth rate compared with C-mesh. Generally, variations of the accuracy of hydrodynamic coefficients are consistent with the increase in the number of cells, and relatively greater improvements can be observed for C D and C M taking log R G N c i0 = 2 as the dividing line. For the NOC less than 3.00 × 10 5 , C-mesh is a better choice for predictions of α < 8.22 • when C-mesh and H-mesh share the same order of magnitude. Considering that near-field node distributions of H-mesh are able to capture more precise flow details around the rudder with denser nearfield node distributions, H-mesh with fine mesh resolutions can be applied with enough computational resources. In this part, to balance the accuracy and efficiency, log R G N c i0 = 2 is selected for following simulations on H-mesh, and corresponded uncertainty estimations are shown in Table 10. Table 9. Comparison of results of H-mesh based on various chord-wise spacings with that of experimental benchmark data adopted from [25] in relative differences with 386 nodes in chord-wise direction for reference setting.   For O-mesh in Table 11, the errors of C L and C D are larger than C and H-mesh because the mesh type fails to capture the wake after the trailing edge. Further, variations of layer-wise growth rates are incapable of improving the situation. Therefore, from the perspective of efficiency, further investigations are not carried out on O-mesh considering its obvious defects.
Various chord-wise spacings are applied for Hybrid-mesh, and the results are in Table 12. Compared with that of structured meshes, the calculation accuracy of C L is as good as C and H-mesh, but that of C D is worse than C and H-mesh while still better than O-mesh. C L tends to be stable when the NOC is more than 2.5 × 10 5 , and C D results are gradually improved with smaller chord-wise spacings. In the following part, a chord-wise spacing of 2.23 × 10 −4 c is selected, and validation for C L is achieved as indicated in Table 13. Table 11. Comparison of results of O-mesh based on various chord-wise spacings to that of experimental benchmark data adopted from [25] in relative differences with 296 nodes in chord-wise direction for reference setting.   Table 12. Comparison of results of Hybrid-mesh based on various chord-wise spacings with that of experimental benchmark data adopted from [25] in relative differences.   No obvious variations can be observed for both C L , C D and C M in Table 14. The insensitivity to layer-wise growth rates for C-mesh is probably caused by good chord-wise mesh resolutions. A rate of 1.15 encounters convergence problems, while results obtained from r 1.10 are stable. Therefore, r = 1.10 is selected.
For H-mesh in Table 15, the accuracy of C L and C D has subtle improvements within 1%, while variability of C M is slightly larger. The trend shows that a slower growth rate of mesh height leads to better results, but, remarkably, too small rates may reduce calculation efficiencies and contribute to mesh quadrilaterals with too large aspect ratios, which are disadvantageous to the convergence of CFD results. In this case, r = 1.05 is well balanced in accuracy and computation time.
For Hybrid-mesh in Table 16, the connection between the structured region and the unstructured region is influenced by layer-wise growth rates. When r = 1.03, relatively large aspect ratios in boundary layer regions increases the difficulty of convergence. Therefore, r = 1.10 is selected as a value of transition of convergent conditions.

Impacts of Far-Field Node Distributions
With selected near-field node distributions, nodes in far-field regions need to be determined to match variations of cell sizes near wall boundaries to present global meshing strategies. For structured meshes, similar to Section 5.3.1, nodes in different directions are changed according to characters of their topological structures, which can be expressed as: where N r i and N w i denote the number of nodes in radial and wake directions. For Hybridmesh, cell sizes in far-field regions are varied to change mesh properties. Table 14. Comparison of results of C-mesh based on various layer-wise growth rates with that of experimental benchmark data adopted from [25] in relative differences.   Table 15. Comparison of results of H-mesh based on various layer-wise growth rates that of experimental benchmark data adopted from [25] in relative differences.   Table 16. Comparison of results of Hybrid-mesh based on various layer-wise growth rates with that of experimental benchmark data adopted from [25] in relative differences.  Number of cells against radial nodes vary over a considerable extent from nearly 3.00 × 10 5 to more than 1.00 × 10 6 for C-mesh and H-mesh, as shown in Tables 17 and 18. However, such variation has few impacts on magnitudes of any hydrodynamic coefficients, indicating that sufficient mesh resolution along the rudder chord has successfully solved the flow in regions with great gradients. In other words, relatively few nodes in radial directions are capable of ensuring the calculation accuracy. Therefore, the reference case already achieved mesh independence for both C-mesh and H-mesh along with radial directions, and the case is still utilized for benchmark in the following studies.
Similar to impacts of nodes in radial directions, variations of C L , C D and C M against changing cell numbers induced by cell sizes in wake directions are very small, as shown in Tables 19 and 20. Compared with the trend of accuracy variations in Section 5.3.1, impacts of node distributions in the near-field on CFD results are stronger than that in the far-field, which is reasonable considering the decaying disturbance of the object on the flow field with increased distances from the wall boundary. As long as mesh resolutions in the near-field are handled carefully, fewer mesh elements are needed in the far-field. Mesh elements in unstructured regions for Hybrid-mesh are changed as shown in Table 21, and the magnitude of results variations are the same as C-mesh and H-mesh. Further decreasing element sizes of unstructured meshes are memory-consuming, and divergent cases occur when the size is below 0.0625 c. Therefore, the body element size of 0.125 c is selected.
Tables 22-26 show uncertainty estimations for C-mesh, H-mesh and Hybrid-mesh with varying far-field node distributions. Most uncertainty quantities in Tables 22-25 are relatively small, corresponding to the fact that further increase in radial or wake directions has little impact on the calculation accuracy. Table 17. Comparison of results of C-mesh based on various radial nodes with that of experimental benchmark data adopted from [25] in relative differences with 150 nodes in radial directions for reference setting.   Table 18. Comparison of results of H-mesh based on various radial nodes to experimental benchmark data adopted from [25] in relative differences with 150 nodes in radial directions for reference setting.   Table 20. Comparison of the results of H-mesh based on various wake nodes with that of experimental benchmark data adopted from [25] in relative differences with 300 nodes in wake directions for reference setting.   Table 21. Comparison of the results of Hybrid-mesh based on various element sizes with that of experimental benchmark data adopted from [25] in relative differences.

Impacts of Mesh Types
Simulations of different mesh types with varying mesh properties are conducted in Sections 5. 1-5.3. Generally, the accuracy of C L is sufficient (within 5% deviations compared with experimental data) for engineering applications except for O-mesh, while discrepancies of C D between CFD and EFD results are relatively large under high angles of attack. Physically, the lift and drag coefficient are nondimensional forms of lift and drag force, respectively, while the former is mainly generated by the pressure difference between the windward and leeward surfaces, the inviscid component, and the latter primarily contributed by the frictional shear stress, the viscous component. The precise prediction of viscous stress relies on a more exact forecasting model in regions of separated flow, along with an advanced understanding of rudder hydrodynamics, the development of the discretization, and numerical methods [38]. In other words, it is inevitable to observe excessive error of C D with the RANS method and limited computational resources. However, selecting the appropriate mesh type can improve the accuracy of both C L and C D to a certain extent in the current framework. Therefore, based on the simulation results obtained above, the impacts of different mesh types on rudder hydrodynamics are summarized in this section. To evaluate the efficiency and accuracy of those mesh types, the calculation accuracy of C D versus changing numbers of cells in the near-field and far-field is shown in Figures 18 and 19, while blue dotted lines connect minimum errors under each α and red lines mark selected mesh settings for each type. Since there is a common trend between impacts of radial and wake distributions on calculation accuracy for C-mesh and H-mesh, only radial node distributions are selected to perform the estimation.
In Figure 18, for the NOC below 3.00 × 10 5 , C-mesh shows higher accuracy than other mesh types. However, H-mesh with better mesh resolution around the rudder achieves better performance under large angles of attack. What's more, it can be foreseen that a further increase in nodes may lead to more accurate results if computational power is not an issue. According to the tendency of Hybrid-mesh, the accuracy of C D remained stable, and smaller chord-wise spacing cannot make up for defects of unstructured mesh in solving rudder flows.
In Figure 19, even though variations of C D of Hybrid-mesh are slightly larger than Cmesh and H-mesh, results obtained from these mesh types settle in stable levels with fixed chord-wise spacings. Although making further efforts in far-field mesh resolutions can still improve the accuracy, the consequent mesh elements with small size may dramatically increase computation time, as well as bring difficulties with convergence. C-mesh and H-mesh with NOC around 5.00 × 10 5 are selected, which are about 2.00 × 10 5 less than Hybrid-mesh in the number of cells and show higher accuracy especially under high rudder angles, proving that structured meshes have higher efficiency in computations than hybrid meshes.  (c) Hybrid-mesh. Figure 18. Calculation deviations of C D compared with experimental benchmark data adopted from [25] varying with numbers of cells, while changing chord-wise spacings in the near-field (blue dotted lines connect minimum errors under each α and red lines mark selected mesh settings for each type).
Based on the principle of balancing efficiency and accuracy, recommended mesh properties for four mesh types investigated in Sections 5.1-5.3 are presented in Table 27, and performance of these settings are compared with that of other codes (Langley Research Center [19]) in Table 28, which coincide with H-mesh best. The results of these CFD codes are achieved using the same block-structured C-mesh in a domain of 500 c around the profile and 500 c behind it, which is much larger than the domain in this paper. Hence, the appearance of differences in Table 27 is reasonable, and relatively good results are achieved by the current method efficiently. Structured meshes of C-type and H-type can achieve accurate C L predictions, while O-mesh fails to capture flow features in wake directions because of the difficulties of actualizing clustering grid points in the wake while maintaining a reasonable spacing at the trailing edge [5]. The concentration degree of mesh elements in regions with a large change of solution gradients for H-mesh is greater than that of C-mesh, contributing to higher accuracy in C D predictions when the NOC is above 3.00 × 10 5 . Some unexpected divergent cases can be observed for C-mesh and H-mesh, and no general rules can be derived to explain why that happens, though it probably relates to the processing mechanism of the solver. In systematic parameter studies, it is recommended to fine-tune mesh properties to obtain convergent results. Hybrid-mesh can be generated more simply with less manual intervention, and calculation results are acceptable but still inferior to C-mesh and H-mesh. However, the high adaptability of unstructured meshes to complex geometries extends the scope of application of Hybrid-mesh, which means that such a meshing strategy can be easily transformed into sophisticated configurations like high-lift rudders. 4 (c) Hybrid-mesh. Figure 19. Calculation deviations of C D compared with experimental benchmark data adopted from [25] varying with numbers of cells, while changing radial node distributions in near-field (blue dotted lines connect minimum errors under each α and red lines mark selected mesh settings for each type).  Table 28. Comparison of mesh types with recommended settings to that of other CFD codes adopted from [19] in relative differences.

Further Validations
To further validate the proposed mesh strategy, more 2D rudder profiles are studied in this section, i.e., NACA 0006, NACA 0009, NACA 0015, and NACA 0018, and the calculated force coefficients are compared with corresponding experimental data. The test conditions are shown in Table 29. 3.26 × 10 6 Jacobs and Sherman [40] By applying H-mesh and Hybrid-mesh with recommended domain size and nodes distributions, obtained C L and C D are compared with EFD results before stall for tested rudder profiles as shown in Figure 20. For profiles with relatively small thickness ratios, such as NACA 0006 and NACA 0009, Hybrid-mesh failed to obtain credible force coefficients. As for H-mesh, in such cases, the mesh type can achieve acceptable C L calculations, although the scope of α which leads to convergent solutions is limited. Nevertheless, from the perspective of structural strength, thin profiles are generally not used for rudder designs, so no further efforts are made to modify current mesh strategies to adapt thin profiles in this paper. For NACA 0015 and NACA 0018, when α < 5 • , some divergence may occur for H-mesh, while Hybrid-mesh shows good convergence. After that, CFD results for H-mesh and Hybrid-mesh agree with EFD ones well, while H-mesh shows better accuracy for α around 15 • , demonstrating that the proposed meshing strategy can be used for other rudder profiles.

Conclusions
This paper describes a comprehensive approach to identifying suitable mesh properties for CFD simulations of airfoil-shaped ship rudders. The goal is to build up an efficient mesh with a sufficient resolution. The suggested mesh properties are summarized in Table 30 for different research objectives. Important conclusions drawn from this work include:

1.
The mesh type affects the efficiency of mesh generation and the accuracy of the CFD solutions. The structured mesh gives more accurate drag predictions with fewer cells compared with that of the hybrid mesh. However, an unstructured mesh is more suitable for complex rudder profiles; it can also provide comparable accurate lift predictions with the structured mesh.

2.
The domain size has to be reasonably large to fully develop the flow. The drag prediction is more influenced by the domain size than the lift prediction. An increase in domain size leads to an improvement in prediction accuracy.

3.
Lift coefficients increase while drag coefficients decrease with an increase of Reynolds numbers in a range of 2 × 10 5 ∼1.89 × 10 7 . The drag coefficient is more sensitive to changes in the Reynolds number than the lift coefficient. A Reynolds number of 6.00 × 10 6 is suggested as a benchmark, above which few changes in lift and drag coefficients were found.

4.
The boundary conditions applied in this study are velocity-inlet, no-slip wall, and pressure-outlet. RANS simulation results with these conditions are in good agreement with the published experimental data [25,39,40].

5.
In all cases, C D prediction accuracy remains poor compared with that of C L . We realize that the above conclusions are based on a case study of 2D ship rudders. For unconventional rudder types, such as flap rudders and fishtail rudders, the conclusions should not be applied directly. For instance, the gap of the flap rudder and the concave part of the fishtail rudder may require an additional number of cells. What's more, attention should be paid to mesh elements in span-wise directions for 3D rudders. Considering larger calculation errors with larger angles of attack, extra efforts are also needed in applications of other CFD methods, such as Large-Eddy Simulation (LES) and Direct Numerical Simulation (DNS), which commonly require a much finer mesh than the RANS method. Therefore, further study is required to determine the most suitable mesh for different rudder types and meshing strategies in 3D configurations.