Numerical Flow Characterization around a Type 209 Submarine Using OpenFOAM †

in the 15th OpenFOAM Workshop (OFW15). Abstract: The safety of underwater operation depends on the accuracy of its speed logs which depends on the location of its probe and the calibration thoroughness. Thus, probes are placed in areas where the ﬂow of water is smooth, continuous, without high velocity gradients, air bubbles, or vortical structures. In the present work, the ﬂow around two different submarines is numerically described in deep-water and near-surface conditions to identify hull zones where probes could be installed. First, the numerical setup of a multiphase solver supplied with OpenFOAM v7 was veriﬁed and validated using the DARPA SUBOFF-5470 submarine at scaled model including the hull and sail conﬁguration at H / D = 5.4 and Fr = 0.466. Later, the grid sensitivity of the resistance was assessed for the full-scale Type 209/1300 submarine at H / D = 0.347 and Fr = 0.194. Free-surface effect on resistance and ﬂow characteristics was evaluated by comparing different operational conditions. Results shows that the bow and near free-surface regions should be avoided due to high ﬂow velocity gradient, pressure ﬂuctuations, and large turbulent vortical structures. Moreover, free-surface effect is stronger close to the bow nose. In conclusion, the probe could be installed in the acceleration region where the local ﬂow velocity is 15% higher than the navigation speed at surface condition. A 4% correction factor should be applied to the probe readings to compensate free-surface effect.


Introduction
Although submarines are the most powerful underwater marine vehicles used to guarantee the maritime sovereignty of nations, their operational safety depends on several mechanical and electronic equipment. For instance, the warship underwater position estimation relies on the accuracy of its speed log, which is affected by the location of its probe on the hull and the calibration procedure to compensate the hydrodynamic nonlinearities on the flow velocity caused by the hull curvature. The accuracy of the measured speed depends on technology used by the log, such as electromagnetic (EM), pitometer, Doppler, impeller, and Global Positioning Satellite (GPS). Although latter is the preferred alternative for modern surface vessels due to its high precision, it is useless for submarines operating completely submerged. One of the viable options for submarines navigating completely submerged are EM logs with flush or protruding probes that estimate the local speed of water by measuring the amplitude of the electrical signal and whose accuracy relays on the probe protruding distance from the hull, flow linearity, turbulence intensity at the probe location, and sensor calibration thoroughness.
The procedure to determinate the probe installation location requires to have an insight into the physics of the water flowing around the full-scale submarine's light hull considering every hydrodynamic variable available in deep-water and near-to-surface conditions. These flow characteristics depend on fluid properties, hull geometry, appendage distribution, vessel motion orientation, and free-surface proximity. Typically, towing tank tests are used to measure the total resistance force, local pressure and velocity, and streamline patterns around underwater vehicles at model scale. Later, full-scale extrapolation is performed by parts: (i) skin friction and pressure components are estimated assuming an axisymmetric body moving through an unbounded fluid, and (ii) wake-making component is estimated considering wave train generation and non-linear effects from free-surface proximity that modifies the pressure distribution around the light hull. For the last step, it is possible to perform full-scale CFD simulations to characterize this free-surface effect on the flow pattern around submarines through detailed analysis of streamline distribution, boundary layer development, vortex core identification, velocity and pressure distributions and turbulence levels. The probe calibration procedure could include free-surface effects using CFD results after a proper verification and validation procedure [1,2] is performed using experimental data from towing tank experiments or sea trials measurements .
One of the experimental databases used for CFD validation was published by the SUBOFF program sponsored by the Defense Advanced Research Projects Agency (DARPA), where a submarine with polynomial geometry was developed and tested including typical appendages at several operational conditions [3,4]. Thereafter, many studies have been performed to understand the flow field around the SUBOFF submarine through numerical prediction of resistance, maneuverability, self-propulsion, both in deep-water and nearsurface conditions. Regarding deep-water analysis, Gross et al. [5] numerically studied the velocity, eddy viscosity, and skin friction at different angles of attack of the bare hull configuration. They confirmed the presence of streamwise vortices on the leeward side at large angles of attack and identified a laminar separation near the submarine stern at low Reynolds number. A thorough investigation of a more complex geometry configuration was performed by Lungu [6], who analyzed the hydrodynamic performance of SUBOFF hull including the sail and four stern appendages using a modified Detached Eddy Simulation (DES) approach. In this work, the streamwise velocity, the pressure distribution and turbulent structures of the nose of the SUBOFF bow and sail confirmed the expected flow physics and showed that the wake at the propeller plane has significant velocity defects, which suggest that special attention is required when modeling the propeller performance.
Regarding free-surface effects, Vasileva and Kyulevcheliev [7] numerically predicted the submarine resistance at different depths for a wide range of velocities and demonstrated that the submarine resistance increases as its free-surface distance decreases due to the wave train generated around submarine. However, this effect decreases at higher speeds in agreement with the different evolution of wave making and viscous resistance coefficients as a function of Froude number. Another alternative used to assess nearsurface conditions was proposed by Gourlay and Dawson who used a Havelock Source Panel Method [8] to take advantage of its low computational requirements. Wave-making resistance was successfully compared with model tests data results except in extreme near-surface condition.
Additionally, several self-propulsive simulations, where the DARPA SUBOFF submarine is fitted with a seven bladed INSEAN E1619 propeller, have been performed to assess the propeller thrust, torque, advance coefficient and wake fraction. Di Felice et al. [9] investigated the performance and fluid dynamics of the INSEAN E1619 propeller running in open water conditions through experiments and simulations. It was suggested that submarine propellers require more cell elements than conventional propellers to capture the weaker flow structures generated by the propeller motion. For this reason, different methods such as body force and sliding mesh method [10] or LES method are used to simulate the flow [11] and thrust identity method [12]. For example, Drogul et al. [13] demonstrate the robustness of the actuator disc theory on the self-propulsion point estimation by coupling open water propeller results with a RANSE solver in a single phase analysis. Here, two methods of self-propulsion were analyzed; (i) where the propeller was modeled as an actuator disc based on a body force method, and (ii) where the propeller behind the hull is included in the simulation. It was shown that the former underpredicts the propulsive efficiency and overpredicts the delivered power. In a complementary study, Delen et al. [14] compared the hydrodynamic performance predictions of DARPA SUBOFF with the previous actuator disk method and empirical method [15]. The CFD resistance predictions are more accurate than those of the empirical method when compared with the experimental data. Finally, Sezen et al. [16] successfully extended the resistance and selfpropulsion analysis to bare hull, sail, and four stern appendages configuration, modeling the propeller action with a Moving Reference Frame (MRF) method. It was found that the hull efficiency is higher in self-propelled case suggesting that this could be due to positive propeller-hull interaction.
Regarding maneuvering analysis, Duman et al. [17] investigated the propeller effects on maneuvering forces, moments, and its derivatives in six oblique towing conditions of the bare hull configuration approximating the propeller effect through an actuator disk. It was found that using this actuator disk model, the propeller does not affect the sway force and yaw moment, but produce a considerable difference in longitudinal forces at small drift angles. Moreover, Chase et al. [18] performed self-propulsion simulations of the fully appendaged DARPA Suboff traveling straight-ahead and turning with several drift angles using RANS and DDES approaches implemented in the CFDShip-Iowa software efficiently coupled with a modified version of a vortex-lattice potential flow propeller code PUF-14 with and a discretized propeller. The most complex case considered was a surfacing maneuver simulation of the DARPA and a discretized E1619 propeller including a sea state of 4. It was showed that the submarine adopts a bow-down attitude as it reaches the surface, that is accompanied with large fluctuations in thrust and torque because of the highly unsteady inflow to the propeller. Furthermore, this configuration was used to analyze a methodology used in surface ships, concluding that this could be applied to compute selfpropulsion of submarines [19]. The comparison of four different turbulence models (RANS, DES, DDES, and no turbulence model) showed that RANS overly dissipates the wake and no turbulence model tip vortices become unphysically unstable. The results confirmed that the approach applied to surface ships is applicable to self-propulsion performance prediction of submarines. However, to the best of the author knowledge, there is no hydrodynamic study of the flow around a different class type of submarine, other than SUBOFF geometry, considering its full-scale length.
Additional CFD studies on the free surface effect on underwater vehicles with different submarines types, other than SUB-OFF, have been performed. For example, the total drag and wake of an axisymmetric UWV model were studied at different depths and speeds using ANSYS CFX [20]. Numerical predictions showed good agreement with the experimental data, where drag coefficient is inversely proportional to Reynolds number for all submergence depths and it is inversely proportional to submergence depth form constant Reynolds number due to the wave-making resistance. In addition, Moonesun and Korol [21] analyzed two torpedo-like models with different L/D ratios using Flow Vision and found that the free surface effects banished for high Froude numbers at an immersion depth of 4.5 D. For regular and irregular waves, a torpedo-like model was analyzed using Flow 3D and Panel method respectively [22,23]. It was recommended an immersion depth of 10% of the incident wavelength to operate safely for AUVs, research submersibles and submarines. Furthermore, a benchmark geometry model, proposed by the CSSRC, was used to assess the hull/propeller interactions in submergence and near surface conditions [24]. Numerical results of thrust, torque, and self-propulsion factors showed a fairly good agreement when compared with experimental data with a small variation of the self-propulsion factors in the near-surface condition. An additional experimental and numerical study on hull/propeller interaction of a realistic modern geometry and a five-blade Wageningen B-series propeller was done by Vali et al. [25]. It was found that when the propeller is operating near the free-surface, the thrust deduction and wake factors decrease and that the thrust coefficient is reduced by 9% when compared to the open water condition. Moreover, Carrica et al. [26] simulated the near-surface selfpropulsion of the generic Joubert BB2 submarine considering two conditions, even keel and in controlled free running considering calm water and regular waves representative of sea state 2 to 7. Here, the top of the sail distance to the calm water surface is used as a reference depth. It was shown that near surface operation causes an increase in the propeller thrust and in waves presented a decay of the wave influence with depth and with decreasing wave amplitude. On the other hand, a compensation through trim tanks is required in controlled free running in response to the vertical forces and pitching moments. Another study was carried out by Zhang et al. [27] and included experimental and numerical simulations of the emergency rising of a Russia's Kilo-class diesel-electric submarine . The effect of course keeping on the roll stability was assessed and it was demonstrated a coupling effect between roll and yaw. It was concluded that the excessive yaw deviation is the cause of emergency rising roll instability at moderate incidence angles.
Among the few available studies on full-scale submarines operating in fully submerged condition, Moonesun et al. [28] analyzed a torpedo-like and an IHSS model through CFD, and a Tango nose model in a towing tank. In all three cases, it was shown that the total resistance coefficient becomes Reynolds number independent after 5 × 10 6 , suggesting that there is a critical Reynolds number that is independent of the submarine geometry. Finally, full-scale numerical simulations of a modern submarine type were done to estimate the hydrodynamic derivatives required to evaluate its maneuvering characteristic [29]. Simulations considering the submarine moving in a straight line and rotating showed that asymmetric forces are generated by a vortex flow developed behind the conning tower that affects diving motion. Hence, all the previously described works focused on global variables estimation using experimental or numerical methods, including hull-propeller interaction without describing the flow structure on the near bow region where speed sensors are usually fitted.
In this work, the flow around a Type 209 submarine is simulated at full-scale using OpenFOAM to identify the suitable hull regions for installing speed logs to guarantee sensor precision. First, the numerical setup of a multiphase solver supplied with OpenFOAM v7 was verified and validated using the DARPA SUBOFF-5470 submarine at scaled model including the hull and sail configuration navigating at H/D = 5.4 and Fr = 0.466. Later, the grid sensitivity analysis of the resistance force is assessed following the Validation and Verification guidelines recommended by the ITTC. Afterward, several fluid characteristics such as: dynamic pressure, velocity on boundary layer limit, skin friction, and vortex cores are used to understand the flow evolution around the submarine light hull at the emerged condition. Later, the influence of the free surface is evaluated using the validated mesh density to identify the flow variations, especially on velocity distribution, including a fully submerged condition. Finally, the linear variation of local flow speed at the probe installation location is assessed to confirm measurement usefulness. It is expected that the methodology used in this work to determine appropriate areas on the hull for probe installation could be generalized to existing or new submarines type designs.

Numerical Method
The air and water flow around the submarine light hull was modeled as a steady-state multiphase flow using OpenFOAM v7 [30]. Following current numerical hydrodynamics practices used for surface vessels, this study adopted a Volume of Fluid (VOF) method, implemented in the interFoam solver, to capture free-surface evolution when the submarine is navigating at near-surface and deep-water depth conditions. Also, the κ − ω SST eddy viscosity model including wall functions were used to capture turbulent structures. The interFoam solver approximates in each cell the Continuity and the Navier-Stokes equations for two immiscible, incompressible, and viscous fluids. An additional transport equation is included to capture free-surface evolution as follows ∂ρ ∂t where ρ is a mixture density defined as ρ = αρ 1 + (1 − α)ρ 2 , being α = 1 if cell is filled with water, and 0 otherwise. Hence, the free-surface interface is located when 0 < α < 1, and the available MULES correction algorithm is used for the free-surface reconstruction guaranteeing its sharpness. In general, equations are integrated in time using the lo-calEuler scheme with a global maximum CFL number of 4 even at the free surface interface. At post-processing, the local flow velocity (V f low ) over the hull and outside boundary layer is estimated from dynamic pressure (p_rgh) and nondimensionalized using submarine navigation velocity (V sub ), as follows: 2.1. DARPA SUBOFF and Type 209/1300 Submarine Geometry As described above, the geometry of the DARPA SUBOFF model is defined by polynomial equations [31] and consists of an axi-symmetrical hull tested in eight different arrangement configurations (AFF) depending on the appendages included [3]. In the present study, only the geometry of the bare hull with the sail is modelled because log probes are typically fitted on the forward half of the hull to avoid stern appendages or propeller interaction. Hence, the experimental resistance force data measured in the DTMB Tow Tank that is used in the validation is estimated from AFF-1 (bare hull), AFF-3 (hull with four stern appendages), and AFF-8 (hull with sail and four appendages).
On the other hand, Type 209 submarine is a diesel-electric warcraft built by HDW of Germany since 1971, with 61 units still operating in 13 countries worldwide. There are five variants depending on the surface displacement, from 1100 to 1500. In the last decade, most of the oldest submarines have being refitted including electronic equipment upgrades. In particular, the 209/1300 submarine analyzed in this study upgraded to two EM speed logs, located at starboard and upper deck, see Figure 1, that must rely on experimental or numerical hydrodynamics studies to determine its light hull location and further calibration curve for near-surface and deep-water navigation. This calibration curve is estimated from sea trials to guarantee the accuracy of the readings reported by the probes. This curve compares the submarine navigation speed and local flow velocity at probe location to compensate non-linearities for hydrodynamic variations produced by the curvature of the hull. In the submarine case, two calibration curves, near-surface and deep-water operation conditions, are used. Typically, the calibration of the deep-water curve uses the near-surface one as the reference navigation speed. However, this heuristic methodology wrongly assumes that flow at deep-water and near-surface conditions are similar, despite experimental evidence proving otherwise, as previously discussed. One alternative to quantify free surface effects at the full-length scale at regions near probe locations is to use an open-access experimental database, such as DARPA SUB-OFF. However, there are significant differences in displacement and length proportions with the 209/1300 class, as shown in Table 1, that prevent a straight forward extrapolation. Hence, the only feasible option is to perform a numerical study at full-scale to characterize flow conditions around the Type 209/1300 submarine light hull. Finally, Table 2 shown the coordinates of four points, P1 to P4, that corresponds to the base and measurement location of fixed protruding probes at starboard (P1, P3) and upper deck (P2, P4), as shown in Figure 1. The distance between each pair of points is 40 cm, corresponding to distance between the probe base and its measurement position. This control points are used to assess the lineal variation of the local flow speed with navigation speed.

Simulation Setup for Submarines
For the DARPA SUBOFF-5470 submarine, the bare hull and sail configuration at model scale was considered, as shown in Figure 2. The hull is Axis-symmetric and is defined by polynomials as described in Groves et al. [31]. The 3D CAD geometry was generated using the Grasshopper plugin in Rhinoceros and the surface was exported to OpenFOAM as an STL file. Additionally, after transforming the geometry to mesh type, the water tightness of the geometry was verified with the software Blender. For the Type 209/1300 submarine, a complete CAD model was generated by combining original blueprints provided by shipowner and measurements taken during a scheduled dry-dock maintenance. However, appendages with dimensions smaller than 10 cm such as: hydroplanes, passive sonars, antennas, and collision bumpers, were not included in the CAD version used in the current work. This simplified CAD version, shown in Figure 3, was used in the numerical simulations to decrease the required computational time per simulation, where neither propeller geometry nor its hull interaction is included. Although the flow around the submarine may be vortex dominated, the dominant vortical structures are generated by the sail [32], that produce asymmetrical flow structures about the central plane that are shed downstream along the hull, hence requiring a fulldomain approach. Nevertheless, there is evidence that when a submarine is moving straight-ahead, it is possible to use the half-domain approach [5,7,10,21]. On the other hand, if the submarine propeller [9,17], self-propulsion problem [12][13][14]18,26], or maneuvering characteristics [27,29] are numerically assessed; then a full-domain approach must be used. The present work focus on the forward hull region below the free surface interface, where the flow is expected to be axisymmetric when fully submerged as reported by Lungu [6], that is slightly perturbed by the occurrence of breaking wave in front of the bow due to vertical, cylindrical bow that is symmetric around centerline when operating at surface depth, as reported by Jasak et al. [33]. Thus, half-domain approach is chosen in this work and decreasing the required time per simulation.
Hence, the computational domain size was defined following ITTC guidelines [34], where the inlet, side, bottom, atmosphere, and outlet boundaries are placed far enough from the submarine to avoid spurious flow interaction, see Figure 4. The boundary conditions are shown in Table 3 for all patches of the computational domain, where the no-slip wall condition is imposed on the submarine boundary. Regarding initial conditions, the flow velocity is initialized with the navigation speed in the inlet, bottom, and internal field.
In addition, turbulent kinetic energy (κ) and its specific rate of dissipation (ω) are set assuming a turbulence intensity level of 1.5% considering wall functions for the hull patch.  The computational grids are generated in stages: (i) an hexahedral background grid using blockMesh with refinement in the region where free-surface interface is expected, (ii) six refinement levels using topoSet tool to half the element size progressively as the elements get closer to the submarine hull, (iii) the surface is inserted with additional (1 or 2) refinement levels if the elements is less than 0.1 D from the submarine hull, (iv) eight layers of prismatic elements parallel to submarine surface are included with a target yPlus, that is estimated considering the laminar boundary layer thickness as a reference length to capture boundary layer evolution efficiently. Figure 5 shows slices of centerline showing one of the generated grids for SUBOFF and S209 submarines. The progressive refinement in regions close to the submarine is evident, aiming to capture high gradients of velocity in the boundary layer and pressure distribution around the submarine's light hull.  Table 4 shows the number of cells used for modelling the computational domain to estimate grid errors and uncertainties. The element size in the inner topoSet region is used as the reference length (∆x) and it is uniform in all three directions to improve the reliability of free surface evolution and turbulence prediction. The Medium grid was generated using the ITTC recommendations and the additional grids were generated using a refinement ratio of rG = √ 2 in both directions. For the Medium grid: (i) DARPA SUBOFF, an element size close to the submarine of 0.05 m (L/∆x = 87.12) and a target yPLus of 50 were used, (ii) Type 209/1300, an element size of 0.20 m (L/∆x = 291.45) and a target yPLus of 300 were used. For the second geometry, the element size is of the same order of magnitude than the probe dimension and it will capture the relevant flow characteristics for speed log accuracy.  Figure 6 shows the distribution of the number of surface layers obtained for each of the considered meshes, despite having requested eight layers, suggesting that most of the hull has the necessary elements to capture the viscous effects within the boundary layer. This difference is due to the mesh quality parameters and curvature of the light hull at the bow and stern sections. Despite the above, the prism surface layer reaches the requested thickness in most of the outer hull, not shown here. Although the yPlus distribution near the surface at the bow, sail, and aft regions has values greater than 300, it is expected that it will not affect the main flow characteristics near the probe location due to its small influence area.

Test Matrix
The DARPA SUBOFF submarine was used for the Verification and Validation considering H/D = 5.4 and Fr = 0.46 condition, where H is the distance from the free-surface to the hull diameter centerline. The numerical prediction is compared to the experimental data provided by Liu and Huang [3]. Additional simulations predictions using the verified grid were compared with the experimental data provided by Neulist [35] for different immersion depth H/D ratios and Froude number, as shown in Table 5: For the Type 209/1300, the navigation conditions shown in Figure 7 were used to assess the flow characteristics, free-surface effects, and potential zones on the light hull of the submarine for probe installation to guarantee measurement accuracy fulfilling the manufacturer requirements. Thus, three different drafts measured at the aft perpendicular were considered without initial trim, namely: (i) surface depth at 5.18 m, (ii) periscope depth at 12 m, and (iii) deep-water at 40 m corresponding to a non-dimensional H/D of 0.33, 1.47, and 6.14 respectively.    Table 6 shows the navigation details for each Type 209/1300 simulation performed in the present study, such as velocity, draft, and grid density considering operational limits reported by the submarine crew. First, simulations 1-3 were used to assess grid convergence independence with the submarine navigating at 9 knots (Fr = 0.194) at surface depth and no trim. Later, simulations 7-9 and 12-14 are completed to confirm grid convergence study navigating at 6 knots (Fr = 0.129) at periscope and deep-water depths respectively. Additional simulations with the submarine at conditions 1, 2, and 3 were used to verify the linear variation of the local speed measured by the probe including a speed of 3 knots (Fr = 0.065).

Results
The numerical results obtained with the settings described in the previous section are summarized here. First, a Verification and Validation procedure is performed using the DARPA SUBOFF submarine geometry at model scale. Next, the grid sensitivity at full-scale was assessed using the total resistance predictions acting on the Type 209/1300 submarine sailing at three navigation conditions: surface depth, periscope depth, and deepwater immersion depth, considering the recommendations described by [1, 33,36,37]. Later, the flow was characterized considering the fine grid by assessing the flow around the external hull, wave train, dynamic pressure, yPlus, and wall Shear Stress distribution, and visualization of vortex cores following [38,39].

Verification and Validation
For this procedure, simulations were performed considering the navigation conditions, H/D = 5.4 and Fr = 0.466, previously described in Table 5. Figure 8 shows: the evolution of: (i) the residual for dynamic pressure, the fraction of water, κ, and ω variables and (ii) viscous and pressure force components acting on the hull, both for the previously described navigation condition using the fine grid with 4,500,437 cells. All residuals smoothly decreased several orders of magnitude, but the dynamic pressure oscillates around 10 5 . Also, both force components achieved convergence after 2000 time steps, which could be related to the absence of free-surface effects. The others grids show the same behaviour.  Figure 9 shows the evolution of the resistance force for all grids. It is evident that convergence is achieved after 2000 iterations. Additionally, all curves are smooth, but Grid 03 that has spurious oscillations. Grid 01 and 02 are close, suggesting that grid convergence have been achieved. Table 7 shows the average values of pressure and viscous force components acting on the whole hull area, and it is estimated as two times the average of the last 1000 iterations for each grid considered in this study. Although the numerical prediction approximates to the experimental value and the error percentage decreased considerably from Grid 05 to 01, it is still at 4.7% for the latter. It is believed that a possible explanation of the current simulations is due to the existing limitations of snappyHexMesh to add enough prism layers in about 1.1% of the hull area, especially on the hull-sail intersection and the aft end of the hull. Another source of error is from the procedure estimation of the experimental data, since the hull-sail configuration is obtained from three different configurations (AFF1 + AFF8 -AFF3) as described in [3]. Further simulations of SUBOFF use the fine grid due to current hardware limitations.   Tables 8 and 9 show the results of the Verification calculations obtained using the Total Resistance coefficient from the numerical results previously shown in Table 7. The convergence ratio R G of 0.172 and 0.605 means that the numerical results for these grid sets (1-3, and 3-5) are in the convergence region and Richardson extrapolation can be used to estimate the numerical force coefficient corresponding to a grid with infinity number of cells. For the grid 1-3 set, the low values of δ G = 0.8% and U GC = 0.6% suggest that numerical results achieved grid convergence and it could be considered as verified. Finally, Figure 10 compares numerical Total Resistance Coefficient prediction to the experimental values provided by Neulist [35] and Liu and Huang [3]. Although the numerical prediction seems to slightly over-predicts the Total drag coefficient for all H/D ratios but for the lowest Froude number Fr = 0.132, simulations are able to capture the non linear variations at intermediate Froude Fr = 0.31 for the near-surface condition H/D = 1.1. Furthermore, the difference at the lowest Froude condition could be explained by limitations on computational resources and should decrease further refining the grid.  Figure 11 shows the evolution of the viscous and pressure force components acting on the submarine light semihull for simulation 03 using CFL = 4 in the fluid domain and free-surface interface taking advantage of the implicit MULES scheme. As usual, the viscous component is related to the shear stresses of both fluids and the pressure component is linked to normal stresses due to the interaction of fluid motion and waves generated around the submarine. Both resistance components reach steady values after 10,000 iterations, with a smoother behaviour for the coarser. Furthermore, the viscous curve converges faster, and it is smoother than the pressure curve, as expected. The force evolution predicted using coarse and intermediate grids shows a similar behavior but with an increased dispersion on the pressure curve, that is suspected to be produced by smaller vortices being captured, and spurious instability produced by high CFL number. For this reason, the smaller CFL number of 4.0 was chosen to avoid the occurrence of numerical instability. In summary, the viscous component represents roughly 30% of the total resistance force, which agrees with the expected predictions of cylindrical bluff bodies moving near a free surface at moderate Froude numbers. The viscous and pressure force predictions shown in Table 10 correspond to average values of the last 1000 iterations for each of the grids considered in Simulation 03, namely, H/D = 0.33 and Fr = 0.194. The standard deviation is used as an indicator of the variation between iterations of the total force, having a maximum value of almost 460 N (2% of average value) of the pressure component for the finer grid, that is similar to the behaviour described in Jasak et al. [33] where a transient simulation of a general cargo carrier shows the resistance force oscillating due to breaking waves in front of the bow. Table 10. Average and standard deviation of force components vs grid for submarine at condition 1 with Domain CFL = 4 and Free surface CFL = 4.  Figure 12 shows the iteration evolution of the total resistance for grids 01 to 05. All curves seem to achieve a quasi-steady state value with oscillations that could be related to the bow breaking wave shown in Figure 13. In this study, convergence condition (R G ) is calculated as 0.18 for the Grid set 2-4 which means that the solution is converged monotonically. For Grid sets 1-3 and 3-5, R G is calculated as −4.366 and −0.524 respectively, which means that the solution has an oscillating convergence. The Grid Convergence Index (GCI) for the set 2-4 is 1.0%, 0.3%, and 3.9% for sets 1-3, 2-4, and 3-5, where the latter has the higher uncertainty. Furthermore, the total resistance is extrapolated from set 2-4 using the Richardson extrapolation procedure [37] to 24,458.84 N for an infinite number of elements. Hence, the fine grid is considered as acceptable considering its balance between the resolution of the flow characteristics and the required computational resources. All remaining simulations will use this fine grid with 9,145,512 elements.  Figure 13 shows the wave train generated around the Type 209/1300 submarine at condition 1, namely, H/D = 0 and Fr = 0.194. On the left-hand side, it is evident that two Kelvin wave patterns are obtained, regardless of a spurious green water on the submarine upper deck produced by the instantaneous acceleration imposed as initial condition. Regarding both wave systems: (i) the system beginning at the submarine's bow with a 0.4 m wave amplitude forward stagnation point, and (ii) the small system beginning at the higher aft foil with a 0.2 m wave amplitude. On the right-hand side, the water fraction distribution profile shows that the first wave crest is almost high enough to wet the upper deck area. However, the next wave trough is modified by a small wave system generated at the upper deck and light hull intersection. The wave profile seems to be canceled after this area as the flow moves downstream, in agreement with the theoretical prediction of the first wave resistance hollow.  Figure 14 shows the dynamic pressure and wall shear stress distribution around the Type 209/1300 submarine at condition 1 and Fr = 0.194. On the right-hand side, it is possible to distinguish a stagnation area at the bow of the forward body section as predicted by a Pitot-tube like device and several pairs of high-low pressure areas produced by wave crest-trough propagation. Moreover, there are stagnation areas at the nose of each aft foil. On the left-hand side, the wall shear stress distribution seems smooth, but in the near-surface bow region at the bow wave trough, where high gradients are observed and the presence of turbulent eddies at the free surface level is suspected.  Figure 15 shows the flow velocity distribution relative to the navigation speed of the submarine including contour lines between 0.5 and 1.3. The stagnation area at the submarine bow described earlier is followed by a wide high-velocity gradient area where, in principle, the speed probe should not be fitted. Moreover, a small high velocity region is revealed at the first wave trough described before that could be related to vortices being shed from this location. Although P1, corresponding to the side probe base location (marked by a grey dot), is at the end of the acceleration region on a wide velocity region where the local flow speed is between 1.1 and 1.2 times higher than the navigation speed. Velocity probe could be installed here, because of appendage restrictions, if the local flow velocity remains inside this range for the whole range of operational speeds. Hence, it is required to include the results from Simulations 04 and 05.   Figure 16 shows the local flow velocity as a function of navigation speed at two locations: P1 and P3, corresponding to the measurement position of a protruding probe that is 40 cm away from P1. Both trend lines for condition 01 can be approximated by a lineal function with a slope of 1.15, which suggests that the flow velocity remains in the 1.1-1.2 range described before and that there is not a significant difference in local flow if measured by a flush or protruding probe when the boundary layer effect is neglected.  Figure 17 shows vortex cores at the submarine bow identified by using Q and λ 2 parameters. Both methods detect a dominant vortex structure at the bow wave and its downstream propagation along the light hull. Small vortex structures with low dynamic pressure (high flow velocity) are detected at the free surface level near the first trough wave, as suggested by the wall shear stress distribution shown above. The latter method is able to detect smaller structures at the free surface and upper deck level. Finally, Figure 18 shows the numerical prediction of drag coefficients for total, viscous, and pressure force components as a function of Froude number, including ITTC'57 curve for the submarine at condition 01 and Total Drag Coefficient for SUBOFF at H/D = 1.1. The total coefficient prediction of the Type 209/1300 submarine has a different behaviour than the SUBOFF, and it is closer to the variation expected for surface vessels. Viscous coefficient data compares well with ITTC'57 curve despite its light hull has fuller form than standard surface vessels. Furthermore, the pressure component, which includes wavemaking and form factor effects, increases with Froude number as expected, and becomes dominant for higher Froude number. The viscous coefficient difference could be related to insufficient elements being used to describe the hull surface for lower Froude numbers.  Figure 19 shows the wave profile at midplane generated by the Type 209/1300 submarine navigating at H/D = 1.47 and Fr = 0.194. On the left-hand side, the sail deforms the free surface above it due to a Venturi effect and produces a short length wave system beginning at the sail end regardless of no free surface piercing. On the right-hand side, the wave elevation demonstrates a weak free-surface interaction with the submarine light hull, where a small wave amplitude system is generated.   Figure 20 shows the flow velocity distribution on light hull relative to the submarine speed including contour lines between 0.5 and 1.3 for H/D = 1.47 and 6.14. In both conditions, the stagnation area at the submarine's bow, sail, and aft foils are followed by a high-velocity gradient where the probe should not be installed. On the other hand, P1-side probe base, marked by a grey dot, is located on a wide velocity region between 1.1-1.2 times the navigation speed and P2-upper deck probe base in on a slower region between 1.0-1.1 times the navigation speed is used when the submarine is completely submerged. Hence, speed probes could be mounted on both locations given the local flow velocity remains this range for the whole range of operational navigation speeds regardless of its difference with the condition 01 H/D = 0.33. Moreover, the calibration procedure used for the upper deck probe is discouraged because it is located on a different velocity range. Finally, the flow velocity distribution is similar for H/D = 1.47 and 6.14 but at the sail produced by its free surface proximity where changes are evident at the sail side. On the left-hand side, Figure 21 shows local flow velocity as a function of navigation speed at two locations: P1 and P2 for H/D = 1.47. Both trend lines can be approximated by linear functions with a slope of 1.12 and 1.03 respectively, suggesting that the local flow velocity around the side probe at P1 remains in the 1.1-1.2 range predicted for the surface depth condition and that it is safe to use this probe when the submarine is operating at near-surface condition. On the other hand, the upper deck probe at P2 is within a different speed range 1.0-1.1, suggesting that the calibration procedure based on the side probe is not appropriate producing an error of 0.25 knots when the submarine is navigating to 9 knots (Fr = 0.194). On the right-hand side, the local flow velocity at P1 is shown as a fraction of the navigation speed. The local flow velocity decreased almost 4% when the submarine move from emerged (H/D = 0.33) to the deep-water condition (H/D = 6.14) slightly increasing the error of the measured velocity by side probe sensor.

Free-Surface Depth Influence on Flow Characteristics
Finally, Figure 22 shows drag components coefficients as a function of Froude number, including the ITTC'57 curve for submarine at H/D = 1.47, and H/D = 6.14 conditions and the Total Drag Coefficient for SUBOFF at H/D = 1.1 and H/D = 3.3. On the righthand side, viscous coefficient data compares well with the ITTC'57 curve, despite the wetted area includes the sail contribution. The pressure component is almost parallel to the viscous one, suggesting that the form factor could be estimated. Total drag coefficients have the same trend as the experimental values for SUBOFF geometry at H/D = 3.3, where some free surface effects are still present. On the left-hand side, viscous data have the same behavior as before, and pressure data seems to increase as the Froude number increases. The Wave-making coefficient could be estimated by deducting the form factor previously estimated.

Conclusions
Results for DARPA SUB-OFF geometry demonstrate that the numerical methodology used in this work is reliable at the engineering level when 4,500,000 cells are used in the simulations. The Verification and Validation procedure was used to estimate the numerical uncertainty of three grid sets combining 5 different densities. Also, strong nonlinear effects due to proximity to free-surface can be captured for Froude number of Fr = 0.3 at H/D = 1.1.
Considering the results and discussion shown in the previous section, it is concluded that the side probe sensor could be installed at P1 location considering that this area presents a continuous flow without recirculating zones or vortex cores for all the navigation conditions considered. However, turbulent structures are concentrated in the near-surface regions. This probe will be located close to the end of the acceleration region at bow, where the local flow speed is still 15% higher than actual submarine speed. This location is recommended to be used when submarine is at surface condition because the flow is considered linear because its velocity remains at the same level for navigation velocities between 0 and 9 knots. Moreover, speed measurements could be useful for near-surface and deep-water conditions if a 4% correction is included to compensate for free-surface effect. On the other hand, the current calibration procedure used for the upper deck probe at P2 using the side probe measurement at P1 as a reference should be further modified to take into account the difference in the flow velocity range of both locations P1 and P2.
Finally, it is recommended to verify the effect of boundary layer thickness on speed measurements of the flush probe type to guarantee the expected accuracy of the flush probe and to generate validation data during sea trials, including at least 8 calibration points for the flush probe type to compensate the manufacturer declared difference in precision (±0.2 knots) between flush and protruding probe types.
This research could follow several interesting paths: (i) to improve results reliability, authors could perform transient simulations, include ramping initial conditions, improve the prism layer generation to cover 100% of the submarine hull, include a wave damping region at the outlet boundary to avoid reflections; (ii) to assess the non-linear free surface effect for Froude numbers between 0.2 and 0.4, to assess the boundary layer thickness effect on speed measurements or assessing the operational trim effect on flow characteristics. Finally, it will be challenging to assess the Type 209 submarine in self propulsion or maneuvering conditions.

Data Availability Statement:
The data presented in this study are available on request from the corresponding author. The data are not publicly available due to its potential impact on submarine operational safety.
Acknowledgments: This research was carried out using the research computing facilities and/or advisory services offered by Scientific Computing Laboratory of the Research Center on Mathematical Modeling: MODEMAT, Escuela Politecnica Nacional-Quito. Special thanks to Mercy Anchundia for her support and the help provided to solve remote access issues. In addition, Authors are deeply thankful with the Submarine Squadron of Ecuador and the submarine crew for their predisposition to share their know-how.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analyses, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Abbreviations
The following abbreviations are used in this manuscript: