Large Eddy Simulation of the Flow around a Generic Submarine under Straight-Ahead and 10 ◦ Yaw Conditions

: Aiming towards a better understanding of the ﬂow ﬁeld around a fully appended Joubert BB2 submarine model, and in order to complement the experimental investigations of the wake of the hydroplanes and sail, large eddy simulation (LES) with the dynamic Smagorinsky model was conducted. Three sets of grids with a maximum grid number of up to 228 million were designed to perform the LES simulation for the Joubert BB2 under 10 ◦ yaw conditions, with a freestream Reynolds number based on the local freestream velocity and a hull length of Re L = 2.2 × 10 7 . Comparisons of the wake of the cruciform appendage were made with experiments to verify the computational accuracy and to examine the inﬂuence of the spatial resolution. A satisfactory result was more representative of the experiments with the improvement in grid spatial resolution. The evolution characteristics of three co-rotating vortices originating from the cruciform appendage under the most reﬁned grid arrangement are further described in detail under straight-ahead and 10 ◦ yaw conditions. The comparison results show that, in the core-ﬂow region, the resultant velocity, vorticity magnitude, and TKE were stronger and the wake was more complicated under 10 ◦ yaw conditions. Tip vortex tracking under 10 ◦ yaw conditions exhibited signiﬁcant three-dimensional characteristics as the wake developed downstream.


Introduction
Submarines will generate vortex structures with various scales and forms that develop downstream during navigation.Generalized appendages, such as sails, hydroplanes, rudders, decks, and flood holes, create coherent vortex structures, like the horseshoe vortex, tip vortex, Karman vortex, discrete vortex, etc., leading to intense flow vibration.The development of these vortex structures also influences the maneuverability and deteriorates the wake of the platform.The vortex oscillation and its induced pressure fluctuation excite the hull, resulting in flow-induced noise, which seriously influences the stealth [1][2][3].The investigation of the flow physics of a submarine can help us to better understand maneuverability limitations and flow-induced noise sources.
A certain amount of the available literature has described the flow around a submarine, within which experimentation plays an important role, especially in early studies.Fu utilized a PIV (particle image velocimetry) system to characterize the flow field around a sting-mounted captive ONR (Office of Naval Research) Body-1 submarine model in a steady turn [4].Jimenez and Ashok utilized a hot-wire system and an SPIV (stereo particle image velocimetry) system to measure the flow field around an axisymmetric DARPA (Defense Advanced Research Projects Agency) SUBOFF (Submarine Technology Program Office) model, respectively, and flow field experimental databases with different Reynolds numbers, pitches, and yaw angles were obtained [5][6][7][8][9][10].Similar research works on the SUBOFF model were conducted by Khan [11] and Shokrallah [12].
Another DSTO (Defense Science and Technology Organization) generic submarine model with a large deck, sail, and an X-form rudder arrangement [13,14], which can provide a useful representation of a conventional submarine, has been widely discussed in the available literature.The flow field measurements of the fully appended model both under straight-ahead conditions and during a 10 • side-slip were conducted by the Defense Science and Technology Group, with data collected using pressure probes, PIV, and flow visualization of the wool-tuft streamers [15][16][17][18][19][20].It provided a benchmark prototype with which to validate and improve numerical simulations of submarine wakes.
Since a new vertical sail and two horizontal hydroplanes, known as "Joubert BB2", were recently designed to improve stability and control characteristics [21,22], the windtunnel experiment of the model with 10• yaw was conducted again, with China-clay visualization and a high-resolution SPIV system [23][24][25].The wake of this cruciform appendage can be made available to assist validation studies.
With the development of computational fluid dynamic (CFD) methods, the numerical simulation of the flow field around a fully appended submarine model can complement experimental investigations and create opportunities to advance the understanding of the flows, as the experiments only obtain limited flow field information.Careful verification and validation studies should be conducted with experimental data.
Four main CFD methods have been developed for predicting the flow field, and direct numerical simulation (DNS) is not practical in engineering predictions because of its enormous demand for computing resources.Reynolds-averaged Navier-Stokes (RANS), DES (detached eddy simulation), and large eddy simulation (LES) are gradually becoming dominant in computing the flow field around underwater vehicles.For the fully appended hull, however, the RANS seems to lack prediction accuracy somewhat, particularly for the second-order statistical moments and the local flow topology, and the DES and its variants lack good representation for the development of wall turbulence.The development of LES can provide very useful insights into the complicated transient nature of the flow, including unsteady wake flow, flow-induced noise, and vibrations [16,[26][27][28].
The LES method was first proposed by Smagorinsky [29], and the applicability of predicting the flow field around submarines has been verified by many researchers.Shi evaluated the predictive capabilities of LES on the flows around a bare hull model by comparing it with experiments, and a good agreement was observed for the pressure and the skin-friction coefficients on the body and the streamwise velocity in the wake [30].Zhou performed LES computations of the SUBOFF submarine model to investigate the wake characteristics of an underwater vehicle with and without a propeller [31].Rosa performed LES simulations of the same model, with a comparison between towed and self-propelled cases, for boundary layer development over both the hull and wake flow [32].Then, LES simulations that focus on the effects of the Reynolds number on the structure of the boundary layer in the stern area, as well as the near wake, were conducted.Results showed that the influence of the junction vortices on the first-and second-order statistics of most of the stern boundary layer is weakly dependent on the Reynolds number [33].The computational ability of LES for the surrounding flow of submarines under maneuvering conditions has been verified.
Zhang studied the numerical prediction approach for the hydrodynamic force and noise of a SUBOFF submarine appended by an AU5-65 propeller using LES and Powell vortex theory [34], and expanded LES to the field of submarine noise prediction and more complex maneuvering conditions, such as submarine propeller interaction and crashback [35][36][37].The results further validated the numerical prediction ability of LES.Kroll performed an LES of the flow over the SUBOFF submarine appended by the DTMB-4381 propeller in forward mode, and crashback, mean flow fields, and propeller load statistics showed good agreement with experiments and previous simulations [38].
In this paper, we focus on the evolution characteristics of three co-rotating vortices originating from the sails and hydroplanes and the development trend of the far wake flow, and provide a comparison of them under straight-ahead and 10• yaw conditions.Three sets of grids with a maximum grid number of up to 228 million were designed to perform LES simulation for the fully appended Joubert BB2 submarine model for 10 • yaw.A comparison of the wake of the cruciform appendage was made with Lee's work [24], to verify the grid convergence and computational accuracy of LES in terms of predicting the flow field around the submarine.Simulations were conducted at a freestream Reynolds number based on the local freestream velocity and the hull length of Re L = 2.2 × 10 7 , which is higher than those in previous work and more representative of a full-scale submarine.Then, the characteristics of the flow around the submarine are described in detail under straight-ahead and 10 • yaw conditions.These computations further elucidate the structure of the flow around the fully appended Joubert BB2 submarine model, and provide an effective complement to the experimental investigations.
The remainder of the paper is organized as follows: Section 1 provides an introduction to the Joubert BB2 submarine model and some information on the SPIV experiments in Lee's work [24].In Section 3, the numerical methodology is explained, including the LES method, computational domain, boundary condition, and the grid arrangements.Section 4 presents the results and discussion, including the validation of the numerical approach and analysis of the evolution of the flow.Section 5 presents the conclusions.

The Joubert BB2 Submarine Model
Figure 1 shows the Joubert BB2 submarine model with a length of L = 3.826 m and a model scale of 1:18.348.The hull was designed as an axisymmetric body of revolution with a length-to-diameter ratio L/2r m = 7.3, and the bow profile is derived from an NACA0018 curve [39], splined to allow the rise in pressure to occur further aft [20].The cross-section of the sail is based on an NACA0022 with a height of 0.080L and a chord length of 0.157L.Two horizontal NACA0015 hydroplanes are assembled on the sail, and the combined span and the root chord were designed to be 0.117L and 0.033L, respectively.The stern control plane consists of four rudders in an X configuration, and the tailing edge is 0.075L away from the end of the hull.Further detailed descriptions can be found in Joubert [13,14], Bettle [21], and Overpelt [22].In this paper, we focus on the evolution characteristics of three co-rotating vortices originating from the sails and hydroplanes and the development trend of the far wake flow, and provide a comparison of them under straight-ahead and 10• yaw conditions.Three sets of grids with a maximum grid number of up to 228 million were designed to perform LES simulation for the fully appended Joubert BB2 submarine model for 10° yaw.A comparison of the wake of the cruciform appendage was made with Lee's work [24], to verify the grid convergence and computational accuracy of LES in terms of predicting the flow field around the submarine.Simulations were conducted at a freestream Reynolds number based on the local freestream velocity and the hull length of ReL = 2.2 × 10 7 , which is higher than those in previous work and more representative of a full-scale submarine.Then, the characteristics of the flow around the submarine are described in detail under straight-ahead and 10° yaw conditions.These computations further elucidate the structure of the flow around the fully appended Joubert BB2 submarine model, and provide an effective complement to the experimental investigations.
The remainder of the paper is organized as follows: Section 1 provides an introduction to the Joubert BB2 submarine model and some information on the SPIV experiments in Lee's work [24].In Section 3, the numerical methodology is explained, including the LES method, computational domain, boundary condition, and the grid arrangements.Section 4 presents the results and discussion, including the validation of the numerical approach and analysis of the evolution of the flow.Section 5 presents the conclusions.

The Joubert BB2 Submarine Model
Figure 1 shows the Joubert BB2 submarine model with a length of L = 3.826 m and a model scale of 1:18.348.The hull was designed as an axisymmetric body of revolution with a length-to-diameter ratio L/2rm = 7.3, and the bow profile is derived from an NACA0018 curve [39], splined to allow the rise in pressure to occur further aft [20].The cross-section of the sail is based on an NACA0022 with a height of 0.080L and a chord length of 0.157L.Two horizontal NACA0015 hydroplanes are assembled on the sail, and the combined span and the root chord were designed to be 0.117L and 0.033L, respectively.The stern control plane consists of four rudders in an X configuration, and the tailing edge is 0.075L away from the end of the hull.Further detailed descriptions can be found in Joubert [13,14], Bettle [21], and Overpelt [22].In Lee's work [24], the experimental results are discussed in a wind axis coordinate system, with the x-axis direction defined as the direction of freestream, as shown in Figure 2, which is also adopted in the paper for a more direct comparison.The velocity field is obtained with an SPIV measurement system at three selected model-length locations, x/L = 0.511, x/L = 0.650, and x/L = 0.815, and the normal direction of these experimental planes is parallel to the free-stream direction.For further details on SPIV setup and measurements, see Lee [24].It should be noted that the model-length Reynolds numbers were chosen to be ReL = 4 × 10 6 and ReL = 8 × 10 6 in those experiments, which are lower than those of this paper.The conclusion has been drawn that similar experimental re- In Lee's work [24], the experimental results are discussed in a wind axis coordinate system, with the x-axis direction defined as the direction of freestream, as shown in Figure 2, which is also adopted in the paper for a more direct comparison.The velocity field is obtained with an SPIV measurement system at three selected model-length locations, x/L = 0.511, x/L = 0.650, and x/L = 0.815, and the normal direction of these experimental planes is parallel to the free-stream direction.For further details on SPIV setup and measurements, see Lee [24].It should be noted that the model-length Reynolds numbers were chosen to be Re L = 4 × 10 6 and Re L = 8 × 10 6 in those experiments, which are lower than those of this paper.The conclusion has been drawn that similar experimental results for the velocity field, turbulence kinetic energy (TKE), cross-stream Reynolds stress, etc., are shown with different Re L values.According to previous research results on the wake field of submarine models, when the Re L increases from 1 × 10 7 to 3.5 × 10 7 , the dimensionless mean velocity changes between 3% and 10%; the change in Re L only influences the amplitude of the fluid variables, not the peak and valley characteristics, and thus the reciprocal validation of computations and experiments is still reliable.

Large Eddy Simulation
The commercial CFD code StarCCM+ by Siemens PLM, based on the finite volume method, was utilized to conduct the LES simulations presented in the paper.The filtered Navier-Stokes equations are as follows: Here,  is the density, v is the filtered velocity, p is the filtered pressure, I is the identity tensor, and b f is the resultant of the body forces.T is the filtered stress tensor due to molecular viscosity, and stress, where S is the strain rate tensor and is computed from the resolved velocity field.The sub-grid scale turbulent viscosity t  must be described by a sub-grid scale model that accounts for the effects of small eddies on the resolved flow, and the dynamic Smagorinsky model proposed by Germano [40] and modified by Lilly [41] is used here.This approach has shown good performance for a variety of complex marine flows of fully appended submarines [31,37,[42][43][44].
As for the discretization of the governing equations, time integration was performed using a semi-implicit, second-order, two-point backward differencing scheme.Convective fluxes were reconstructed using multi-dimensional, cell-limited linear interpolation, whereas diffusive fluxes were reconstructed using a combination of central-difference approximations and gradient face interpolation to minimize the non-orthogonality error.The coupling between the velocity and pressure was achieved by means of the classic PISO (Pressure Implicit with Splitting of Operators) algorithm, and the algebraic multigrid method was employed to accelerate the solution convergence.

Large Eddy Simulation
The commercial CFD code StarCCM+ by Siemens PLM, based on the finite volume method, was utilized to conduct the LES simulations presented in the paper.The filtered Navier-Stokes equations are as follows: Here, ρ is the density, v is the filtered velocity, p is the filtered pressure, I is the identity tensor, and f b is the resultant of the body forces.T is the filtered stress tensor due to molecular viscosity, and T SGS = 2µ t S − 2 3 (µ t ∇ • v)I is the sub-grid-scale stress, where S is the strain rate tensor and is computed from the resolved velocity field.The sub-grid scale turbulent viscosity µ t must be described by a sub-grid scale model that accounts for the effects of small eddies on the resolved flow, and the dynamic Smagorinsky model proposed by Germano [40] and modified by Lilly [41] is used here.This approach has shown good performance for a variety of complex marine flows of fully appended submarines [31,37,[42][43][44].
As for the discretization of the governing equations, time integration was performed using a semi-implicit, second-order, two-point backward differencing scheme.Convective fluxes were reconstructed using multi-dimensional, cell-limited linear interpolation, whereas diffusive fluxes were reconstructed using a combination of central-difference approximations and gradient face interpolation to minimize the non-orthogonality error.The coupling between the velocity and pressure was achieved by means of the classic PISO (Pressure Implicit with Splitting of Operators) algorithm, and the algebraic multigrid method was employed to accelerate the solution convergence.

Computational Domain and Boundary Condition
Figure 3 shows the cylindrical computational domain; the domain has a length of 4L and a radius of L, and extends from L upstream of the front of the hull to 2L downstream of the stern of the hull, to model a fully developed wake.Free-stream boundary conditions are imposed at the inflow and radial boundaries, and a pressure-outlet boundary condition is imposed at the outflow boundary.The no-slip wall treatment is used for shear stress specification.

Computational Domain and Boundary Condition
Figure 3 shows the cylindrical computational domain; the domain has a length of 4L and a radius of L, and extends from L upstream of the front of the hull to 2L downstream of the stern of the hull, to model a fully developed wake.Free-stream boundary conditions are imposed at the inflow and radial boundaries, and a pressure-outlet boundary condition is imposed at the outflow boundary.The no-slip wall treatment is used for shear stress specification.

Computation Mesh
In this subsection, three sets of grids are designed for the convergence study.The simulation domain is discretized using the unstructured hex-dominant grid (trimmer mesher, Octree grid), and a prism layer mesher is used for the generation of a boundary layer mesh with the dimensionless wall distance y + ≈ 1 for all grid sets.The refinement ratio of grids in three directions is 2 , as recommended by the International Towing Tank Conference [45].The corresponding grid numbers are 3.30 × 10 7 , 8.62 × 10 7 , and 2.28 × 10 8 , respectively, and the three grid codes are G1 to G3, respectively.All numerical simulations are carried out by parallel processing in CSSRC (China Ship Scientific Research Centre) with 50 nodes (2400 processors).The time needed to finish a simulation with 228 million cells is about 30 days.
Figure 4 shows a diagram of the basic grid arrangement and the intuitive comparison of grid spatial resolution around the hydroplanes and sail under different grid sets.For capturing the flow field more accurately, a local volumetric control block was established around the submarine to control the surrounding grid size.Care was taken to resolve the initial vortex formation and roll-up of the free shear layer, and to avoid the rapid dissipation of the vortex structure in the wake of the sail, hydroplanes, and X-rudder; local volumetric control blocks were also established around these appendages for more refined volume mesh control.There is an angle of approximately 5.5 degrees between the blocks and the longitudinal axis of the hull in calculations under 10° yaw conditions.

Computation Mesh
In this subsection, three sets of grids are designed for the convergence study.The simulation domain is discretized using the unstructured hex-dominant grid (trimmer mesher, Octree grid), and a prism layer mesher is used for the generation of a boundary layer mesh with the dimensionless wall distance y + ≈ 1 for all grid sets.The refinement ratio of grids in three directions is √ 2, as recommended by the International Towing Tank Conference [45].The corresponding grid numbers are 3.30 × 10 7 , 8.62 × 10 7 , and 2.28 × 10 8 , respectively, and the three grid codes are G1 to G3, respectively.All numerical simulations are carried out by parallel processing in CSSRC (China Ship Scientific Research Centre) with 50 nodes (2400 processors).The time needed to finish a simulation with 228 million cells is about 30 days.
Figure 4 shows a diagram of the basic grid arrangement and the intuitive comparison of grid spatial resolution around the hydroplanes and sail under different grid sets.For capturing the flow field more accurately, a local volumetric control block was established around the submarine to control the surrounding grid size.Care was taken to resolve the initial vortex formation and roll-up of the free shear layer, and to avoid the rapid dissipation of the vortex structure in the wake of the sail, hydroplanes, and X-rudder; local volumetric control blocks were also established around these appendages for more refined volume mesh control.There is an angle of approximately 5.5 degrees between the blocks and the longitudinal axis of the hull in calculations under 10 • yaw conditions.

Results and Discussion
In this section, we initially present results for the qualitative and quantitative comparison between the numerical results and experimental results for 10° yaw with three sets of grids, and then under the most refined grid arrangement.The detailed numerical

Results and Discussion
In this section, we initially present results for the qualitative and quantitative comparison between the numerical results and experimental results for 10 • yaw with three sets of grids, and then under the most refined grid arrangement.The detailed numerical investigation of the flow around the submarine was conducted under straight-ahead and 10 • yaw conditions.are denoted by Exp. 1 and Exp. 2, respectively.It can be observed that the sail-tip vortex formed by the rolling up of the wake can be presented by numerical calculations under all sets of grids.The region of the core flow, defined as a region of vortex flow from its center to its radial location of maximum swirl, is quite similar to the experimental results, while the range of the core flow and low-velocity region presents minor differences under different sets of grids.For the vertical component of velocity, in particular, as the number of grids increases, the distribution of high/low-velocity regions becomes more concentrated and obvious, and the numerical results become closer to the experiments.In addition, the numerical dissipation decreases as the grid number increases, and the capture of the cross-stream Reynolds stress seems to be more refined.Figure 8 shows a comparison of the mean resultant velocity figures < U xyz > /U ∞ , as functions of radial distance r y /L from the vortex center, in the horizontal profiles through the sail-tip vortex.The mean resultant velocity under different grid sets exhibits a considerable difference, especially for the region of the core flow, which presents a corresponding increase with the grid number's continuous growth.The core-flow velocity obtained from numerical calculations exhibits a significant difference at different axial distances from the bow, but is quite close for the experiments.The core flow obtained in the first set of grids (G1) has a mean velocity < U xyz > /U ∞ ≈ 1.12 and < U xyz > /U ∞ ≈ 0.91 at x/L = 0.511 and x/L = 0.815, respectively, with an axial descent rate of 18.8%.It corresponds to < U xyz > /U ∞ ≈ 1.26 and < U xyz > /U ∞ ≈ 1.11, respectively, in G3, with an axial descent rate of 11.9%.This indicates that numerical attenuation inevitably exists in the LES simulations with the dynamic Smagorinsky model of the core flow, which is not advantageous compared to experiments, and can be partly eliminated by the improvement in the spatial resolution of the grid.  .It can be seen that the numerical results become gradually closer to the experimental values as the grid number increases.For the prediction of the peak and valley values, especially, the results obtained in G3 are quite close to those of Exp. 2, with a relatively higher Reynolds number.Overall, the relative errors of the peak and valley values in G3 are 10.6% and 4.3%, respectively, which becomes 13.2% and 14.0% in G2, and 21.6% and 21.1% in G1.It can be concluded that, with the refinement of the grids, the extreme values can be more accurately captured in numerical simulations due to the improvement in spatial resolution.Figure 9 shows a comparison of the vertical component of velocity < U z > /U ∞ .It can be seen that the numerical results become gradually closer to the experimental values as the grid number increases.For the prediction of the peak and valley values, especially, the results obtained in G3 are quite close to those of Exp. 2, with a relatively higher Reynolds number.Overall, the relative errors of the peak and valley values in G3 are 10.6% and 4.3%, respectively, which becomes 13.2% and 14.0% in G2, and 21.6% and 21.1% in G1.It can be concluded that, with the refinement of the grids, the extreme values can be more accurately captured in numerical simulations due to the improvement in spatial resolution.Figure 11 summarizes the comparison of the centerline of the sail wake, which is defined by tracing the boundary <uxuy > = 0 between the positive and the negative Reynolds stresses.Overall, the sail-wake centerline shifts very slightly windward by no more than 0.01L, and by refining the grids from G1 to G3, it shifts leeward and closer to the experiments.The centers of the vortices on the upper hull are listed in Table 1, and their relative errors to Exp. 2 are listed in Table 2.In the near sail-wake region (x/L = 0.511), the numerical simulations for all grid sets can define the centers of the vortices well, but as they evolve downstream, the advantages of refining the grids gradually become reflected.For the case of G3, the maximum horizontal and vertical relative errors of the vertex centers at x/L = 0.815 are 7.2% and 1.5%, which is quite satisfactory considering the actual engineering prediction needs.Figure 11 summarizes the comparison of the centerline of the sail wake, which is defined by tracing the boundary <u x u y > = 0 between the positive and the negative Reynolds stresses.Overall, the sail-wake centerline shifts very slightly windward by no more than 0.01L, and by refining the grids from G1 to G3, it shifts leeward and closer to the experiments.The centers of the vortices on the upper hull are listed in Table 1, and their relative errors to Exp. 2 are listed in Table 2.In the near sail-wake region (x/L = 0.511), the numerical simulations for all grid sets can define the centers of the vortices well, but as they evolve downstream, the advantages of refining the grids gradually become reflected.For the case of G3, the maximum horizontal and vertical relative errors of the vertex centers at x/L = 0.815 are 7.2% and 1.5%, which is quite satisfactory considering the actual engineering prediction needs.

Analysis of the Evolution of the Flow
An overall conclusion can be drawn that the results of the LES simulations are more representative of the experiments as the grid spatial resolution is improved, and then under the most refined grid arrangement, the evolution of the flow under straight-ahead and 10 • yaw conditions is further analyzed.
Figure 12 presents an overall view of the flow past Joubert BB2 under straight-ahead conditions in terms of the second invariant of the velocity gradient Q, colored by the mean resultant velocity < U xyz > /U ∞ .At the immediate front of the junctions of the sail root and the deck, the flow rolls up into a horseshoe vortex system surrounding the sail, the legs of which develop downstream following the deck.Because of the adverse pressure gradient at the trailing edge of the sail, the flow gradually develops into turbulence, and the side vortices are formed and interact with the horseshoe vortex at the root.The side vortices over the sail cap are transported along the sail edge and then merge with the sail-tip vortices and dissipate rapidly.The flow over the outer edge of the hydroplanes induces a pair of hydroplane-tip vortices with opposite circulation, which develop downstream independently, then dissipates and disappears at a distance of approximately one submarine length from the stern.Further, horseshoe vortex, tip vortex, and wake vortex systems can be observed around the X-rudders.Between the two upper rudders, they interact with the vortices sweeping down and sideways over the end of the deck, complicating the flow into the propeller disk, which is unstable and the main source of propeller hydrodynamic noise.Downstream, far away from the hull, all the tip vortices dissipate and really only the wake vortices, after complicated interaction, dynamically evolve, with the energy gradually weakening.
into the wake between the upper and lower rudders.The wake of the submarine becomes quite complicated, and the flow behind the stern is dominated by the mixing of various component vortex systems, including the tilted horseshoe vortex system, the upper and lower hull vortices, the tip vortices, and the wake of the sail, hydroplanes, X-rudders, and hull.Similarly, Figure 13 presents various vortex systems past Joubert BB2 under 10 • yaw conditions, from perspectives of oblique, top, and side views.It can be clearly seen that the vortex systems are more complicated.The side vortices on the leeward side of the sail occur further forward, and the sail-tip vortices are relatively strong enough to develop far downstream.The same is true for the hydroplane-tip vortices, but they are not clearly observed downstream because of strong interactions with the sail wake.Figure 14 shows the development of the trajectories of the cores of tip vortices originating from the cruciform appendage, including the port and starboard hydroplanes and the vertical sail, under 10 • yaw conditions.Referring to Akkermans's work [46], tip vortex tracking under 10 • yaw conditions exhibits significant three-dimensional characteristics as the wake develops downstream.The clockwise rotating sail-tip vortices maintain an axial angle of approximately 8 degrees with the hull and develop downstream and leeward, and are almost stable vertically after experiencing a brief downwash immediately behind the sail.The position of the port hydroplane-tip vortices fluctuates widely, the horizontal and vertical coordinates of which reach their maximum values approximately at x/L = 1.1 and x/L = 0.7, respectively, and then experience a sharp drop.Overall, the port hydroplanetip vortices develop and revolve around the sail-tip vortices.The development of the starboard hydroplane-tip vortices is relatively stable, and their core keeps moving towards the leeward side, with the vertical position gradually rising away from the hull after passing through a valley at approximately x/L = 1.1, due to the repulsive interaction of the hull wake.
Under 10 • yaw conditions, another obvious feature that distinguishes the straightahead conditions is the flow separation on the leeward side of the middle hull.The upper and lower vortex system can be clearly seen, and the former eventually interacts with the horseshoe vortex system originating from the sail, while the latter merges into the wake between the upper and lower rudders.The wake of the submarine becomes quite complicated, and the flow behind the stern is dominated by the mixing of various component vortex systems, including the tilted horseshoe vortex system, the upper and lower hull vortices, the tip vortices, and the wake of the sail, hydroplanes, X-rudders, and hull.Figure 15 shows the evolution of the mean vorticity magnitude < ω xyz > r m /U ∞ along the submarine axial, under straight-ahead and 10 • yaw conditions, where r m = L/(2 × 7.3).The momentum and energy are transported by the development of vortex systems, where the vortex systems generated by the sail, hydroplanes, deck, hull, and X-rudder pass present an increase in vorticity, especially for the tip vortices, horseshoe vortices, wake vortices, and generated hull side vortices, possibly.It can be seen that the instabilities of the horseshoe-vortex system and its interaction with the hull boundary layer cause the legs to break up and develop connected vortex loops, which results in the transport of momentum across the hull and influences the distribution of the vorticity, as mentioned by Fureby [20].As the vortex structure gradually dissipates downstream, the vorticity gradually decreases, with the sail-tip vortices and hydroplane-tip vortices being particularly prominent.For the case of straight-ahead conditions, the obvious vorticity near the sail induced by the tip vortices of the cruciform appendage quickly decays, quite significantly far away from the stern for the case with 10 • yaw.Increased vorticity is also found around and behind the hull under 10 • yaw conditions because of considerable interaction between the flow and the hull.Therefore, the flow-induced noise of submarines under maneuvering conditions has always been a research highlight in the international hydrodynamics field.As the vortex systems develop from a concentrated distribution in the near wake region to a dispersed mode in the far field, in both cases, the vorticity gradually weakens while the coverage range increases, due to energy conservation.Figure 16 shows the evolution of the turbulence kinetic energy 2 / kU  along the submarine axial, under straight-ahead and 10° yaw conditions.The evolution of TKE is quite similar to that of the mean vorticity magnitude, where the vorticity is concentrated, the momentum is intense, and the TKE is significant.For the case with 10° yaw, it can be obviously seen that the TKE induced by the sail-tip vortices is stronger than that induced by the hydroplane-tip, with the same pattern as the vorticity followed.In addition, the TKE in the wake is strongly influenced by the X-rudder under both conditions, which exacerbates the velocity fluctuations of the flow and induces additional propeller noise; thus, the hydrodynamic design of the X-rudder is also a research highlight.Figure 16 shows the evolution of the turbulence kinetic energy k/U 2 ∞ along the submarine axial, under straight-ahead and 10 • yaw conditions.The evolution of TKE is quite similar to that of the mean vorticity magnitude, where the vorticity is concentrated, the momentum is intense, and the TKE is significant.For the case with 10 • yaw, it can be obviously seen that the TKE induced by the sail-tip vortices is stronger than that induced by the hydroplane-tip, with the same pattern as the vorticity followed.In addition, the TKE in the wake is strongly influenced by the X-rudder under both conditions, which exacerbates the velocity fluctuations of the flow and induces additional propeller noise; thus, the hydrodynamic design of the X-rudder is also a research highlight.In view of the fact that the momentum and energy transported by the sail-tip vortices under 10° yaw conditions were predominant, the flow characteristics of the sail-tip vortices are further studied below.Figure 20   In view of the fact that the momentum and energy transported by the sail-tip vortices under 10° yaw conditions were predominant, the flow characteristics of the sail-tip vortices are further studied below.Figure 20  In view of the fact that the momentum and energy transported by the sail-tip vortices under 10 • yaw conditions were predominant, the flow characteristics of the sail-tip vortices are further studied below.Figure 20 provides a comparison of the velocity for the sail-tip vortex under different longitudinal locations, including mean resultant velocity < U xyz > /U ∞ and the three-dimensional component of velocity < U x > /U ∞ , < U y > /U ∞ , and < U z > /U ∞ .It can be seen intuitively that, as the wake of the sail-tip vortices develops downstream, the mean resultant velocity, streamwise velocity, and horizontal velocity show a gradually decreasing trend, and the fluctuation of the vertical velocity between peak and valley values gradually weakens.In the near-wake region, where x/L ≤ 0.807, the core flow exhibits a high-velocity characteristic, while in the far-wake region, the velocity is smaller than the freestream velocity. .It can be seen intuitively that, as the wake of the sail-tip vortices develops downstream, the mean resultant velocity, streamwise velocity, and horizontal velocity show a gradually decreasing trend, and the fluctuation of the vertical velocity between peak and valley values gradually weakens.In the near-wake region, where x/L ≤ 0.807, the core flow exhibits a high-velocity characteristic, while in the far-wake region, the velocity is smaller than the freestream velocity.changes with the evolution of the flow.In the near wake region, where x/L ≤ 0.807, the peaks are located on the windward side and the valleys are located on the leeward side, while in the far wake region, the positions of the two are exactly opposite.Figure 21 shows a comparison of the mean vorticity under different streamwise locations.The dominant component of the mean vorticity magnitude < ω xyz > r m /U ∞ is the streamwise mean vorticity < ω x > r m /U ∞ , which indicates that the sail-tip vortex rotates rather faster around the x-axis compared to the other two.The core-flow vorticity weakens rapidly as the wake develops downstream, while the decay rate gradually slows down.An interesting phenomenon shows that the distribution of the peaks and valleys of the vertical mean vorticity < ω z > r m /U ∞ changes with the evolution of the flow.In the near wake region, where x/L ≤ 0.807, the peaks are located on the windward side and the valleys are located on the leeward side, while in the far wake region, the positions of the two are exactly opposite.
Figure 22 shows a comparison of the turbulence kinetic energy under different streamwise locations.The roll-up from the sail is accompanied by a downwash of the sail-tip vortex (see Figure 14a), and the normal stress in the vertical component <u x u x > is the strongest contribution to the TKE, as mentioned by Lee [24].Because of the significant cross-stream, the normal stress in the horizontal component accounts for the second strongest contribution, and the normal stress in the streamwise component is the weakest.The strongest TKE does not occur immediately behind the sail, but approximately in the range of x/L equal to 0.6 to 1.0, where the downwash of the sail-tip vortex is quite intense.Figure 22 shows a comparison of the turbulence kinetic energy under different streamwise locations.The roll-up from the sail is accompanied by a downwash of the sail-tip vortex (see Figure 14a), and the normal stress in the vertical component <uxux> is the strongest contribution to the TKE, as mentioned by Lee [24].Because of the significant cross-stream, the normal stress in the horizontal component accounts for the second strongest contribution, and the normal stress in the streamwise component is the weakest.The strongest TKE does not occur immediately behind the sail, but approximately in the range of x/L equal to 0.6 to 1.0, where the downwash of the sail-tip vortex is quite intense.

Conclusions
In the paper, large eddy simulation with the dynamic Smagorinsky model was conducted to investigate the flow field around a fully appended Joubert BB2 submarine model under straight-ahead and 10° yaw conditions.Conclusions acquired from the analysis of computed results can be summarized as follows.
(1) Three sets of grids under 10° yaw conditions were designed to examine the grid convergence and computational accuracy.As the grid number increased, the numerical dissipation decreased, and the capture of the vertical component of velocity, the cross-stream Reynolds stress, seemed to be more refined, especially for the extreme values, which were more representative of the experiments, with the relative error of the peak value in the region of core flow being relatively smaller.The numerical attenuation of the mean resultant velocity of the core flow inevitably existed in the LES simulations with the dynamic Smagorinsky model, which could be eliminated through the improvement in the grid's spatial resolution.In the near sail-wake region, the numerical simulations for all grid sets could define the centers of the vortices well, but as they evolve downstream, the advantages of refining the grids are gradually reflected.Overall, through qualitative and quantitative comparison with experiments under 10° yaw conditions, the computational accuracy was

Conclusions
In the paper, large eddy simulation with the dynamic Smagorinsky model was conducted to investigate the flow field around a fully appended Joubert BB2 submarine model under straight-ahead and 10 • yaw conditions.Conclusions acquired from the analysis of computed results can be summarized as follows.
(1) Three sets of grids under 10 • yaw conditions were designed to examine the grid convergence and computational accuracy.As the grid number increased, the numerical dissipation decreased, and the capture of the vertical component of velocity, the cross-stream Reynolds stress, seemed to be more refined, especially for the extreme values, which were more representative of the experiments, with the relative error of the peak value in the region of core flow being relatively smaller.The numerical attenuation of the mean resultant velocity of the core flow inevitably existed in the LES simulations with the dynamic Smagorinsky model, which could be eliminated through the improvement in the grid's spatial resolution.In the near sail-wake region, the numerical simulations for all grid sets could define the centers of the vortices well, but as they evolve downstream, the advantages of refining the grids are gradually reflected.Overall, through qualitative and quantitative comparison with experiments under 10 • yaw conditions, the computational accuracy was verified and results shown are more representative of the experiments with the improvement in grid spatial resolution.(2) A comparison of the evolution of the flow under straight-ahead and 10 • yaw conditions shows that in the core-flow region, the resultant velocity, vorticity magnitude, and TKE under straight-ahead conditions were somewhat smaller than those under 10 • yaw conditions.The side-vortices on the leeward side of the sail occurred further forward, and the sail-tip and hydroplane-tip vortices were strong enough, relatively, to develop far downstream under 10 • yaw conditions.Another obvious feature that distinguishes the straight-ahead conditions is the flow separation on the leeward side of the middle hull; the wake of the submarine becomes quite complicated, and the flow behind the stern is dominated by the mixing of various component vortex systems, including the tilted horseshoe-vortex system, the upper and lower hull vortices, the tip vortices, and the wake of the sail, hydroplanes, X-rudders, and the hull under 10 • yaw conditions.But downstream, far away from the hull under straight-ahead conditions, all the tip vortices dissipated, and only the wake vortices, after complicated interaction, dynamically evolved, with the energy gradually weakening.
(3) The tip vortex tracking under 10 • yaw conditions exhibited significant three-dimensional characteristics compared to those under straight-ahead conditions.Under10 • yaw conditions, sail-tip vortex tracking maintained an axial angle of approximately 8 degrees with the hull, and was almost stable vertically after experiencing a downwash immediately behind the sail.The port hydroplane-tip vortices developed and spiraled around the sail-tip vortices, while the core of the starboard hydroplane-tip vortices kept moving towards the leeward side, with the vertical position gradually rising away from the hull after passing through a valley at approximately x/L = 1.1, due to the repulsive interaction of the hull wake.(4) The resultant velocity, vorticity magnitude, and TKE showed a gradually decreasing trend as the wake of the cruciform appendage developed downstream.Under 10 • yaw conditions, the core-flow exhibited a high-velocity characteristic, in which the peaks of mean vorticity magnitude were located on the windward side in the near wake region, while in the far wake region, the velocity was smaller than the freestream velocity and the valleys of mean vorticity magnitude were located on the windward side.The strongest TKE did not occur immediately behind the sail, but approximately in the range of x/L equal to 0.6 to 1.0, where the downwash of the sail-tip vortex was quite intense.

J
. Mar. Sci.Eng.2023, 11, x FOR PEER REVIEW 4 of 25 sults for the velocity field, turbulence kinetic energy (TKE), cross-stream Reynolds stress, etc., are shown with different ReL values.According to previous research results on the wake field of submarine models, when the ReL increases from 1 × 10 7 to 3.5 × 10 7 , the dimensionless mean velocity changes between 3% and 10%; the change in ReL only influences the amplitude of the fluid variables, not the peak and valley characteristics, and thus the reciprocal validation of computations and experiments is still reliable.

Figure 2 .
Figure 2. The schematic diagram of SPIV measurement and the wind axis coordinate system.

Figure 2 .
Figure 2. The schematic diagram of SPIV measurement and the wind axis coordinate system.

Figures 5 -
Figures 5-7 provide a direct comparison of the wake of the cruciform appendage, including the mean resultant velocity < U xyz > /U ∞ , the vertical component of velocity < U z > /U ∞ , and the cross-stream Reynolds stress < u y u z > /U 2 ∞ , where U ∞ denotes the freestream velocity.The experimental results with Re L = 4 × 10 6 and Re L = 8 × 10 6are denoted by Exp. 1 and Exp. 2, respectively.It can be observed that the sail-tip vortex formed by the rolling up of the wake can be presented by numerical calculations under all sets of grids.The region of the core flow, defined as a region of vortex flow from its center to its radial location of maximum swirl, is quite similar to the experimental results, while the range of the core flow and low-velocity region presents minor differences under different sets of grids.For the vertical component of velocity, in particular, as the number of grids increases, the distribution of high/low-velocity regions becomes more concentrated and obvious, and the numerical results become closer to the experiments.In addition, the numerical dissipation decreases as the grid number increases, and the capture of the cross-stream Reynolds stress seems to be more refined.

Figure 5 .
Figure 5. Qualitative comparison of the mean resultant velocity.Figure 5. Qualitative comparison of the mean resultant velocity.

Figure 5 .
Figure 5. Qualitative comparison of the mean resultant velocity.Figure 5. Qualitative comparison of the mean resultant velocity.

Figure 5 .
Figure 5. Qualitative comparison of the mean resultant velocity.

Figure 6 .
Figure 6.Qualitative comparison of the vertical component of velocity.Figure 6. Qualitative comparison of the vertical component of velocity.

Figure 6 .
Figure 6.Qualitative comparison of the vertical component of velocity.Figure 6. Qualitative comparison of the vertical component of velocity.

Figure 7 .
Figure 7. Qualitative comparison of t velocity cross-stream Reynolds stress.

Figure 8 1 Figure 7 .
Figure8shows a comparison of the mean resultant velocity figures

Figure 9
Figure 9 shows a comparison of the vertical component of velocity

Figure 10 Figure 9 .
Figure 10 provides a plot of the comparison of the turbulence kinetic energy

Figure 10
Figure10provides a plot of the comparison of the turbulence kinetic energy k/U 2 ∞ for the sail-tip vortex.The distribution of TKE along the radial distance r y /L for all grid sets is consistent with the experimental results, and the simulations of G3 are more representative of the experiments, with the relative error of the peak value in the region of core flow being relatively smaller.

Figure 12 .
Figure 12.Development of vortex systems past Joubert BB2 under straight-ahead conditions.

Figure 12 .
Figure 12.Development of vortex systems past Joubert BB2 under straight-ahead conditions.

Figure 12 .
Figure 12.Development of vortex systems past Joubert BB2 under straight-ahead conditions.

Figure 15
Figure15shows the evolution of the mean vorticity magnitude

Figures 17 -, and the turbulence kinetic energy 2 /Figure 16 .
Figures 17-19 show a comparison of the mean resultant velocity

Figures 17 -Figure 16 .
show a comparison of the mean resultant velocity < U xyz > /U ∞ , the mean vorticity magnitude < ω xyz > r m /U ∞ , and the turbulence kinetic energy k/U 2 ∞ , respectively, for the model-length locations x/L = 0.484, under straight-ahead and 10 • yaw conditions.A pair of sail-tip vortices with opposite circulation can be found in straight-ahead conditions, and in the core-flow region, the mean resultant velocity, the mean vorticity magnitude, and the turbulence kinetic energy were considerably smaller than those under 10 • yaw conditions.

Figures 17 - 2 /Figure 17 .
Figures 17-19 show a comparison of the mean resultant velocity

Figure 20 .
Figure 20.(a) The mean resultant velocity; (b) The mean streamwise velocity; (c) The mean horizontal velocity; (d) The mean vertical velocity.Comparison of the mean velocity under different streamwise locations.

Figure 21
Figure 21 shows a comparison of the mean vorticity under different streamwise locations.The dominant component of the mean vorticity magnitude

Figure 20 .
Figure 20.(a) The mean resultant velocity; (b) The mean streamwise velocity; (c) The mean horizontal velocity; (d) The mean vertical velocity.Comparison of the mean velocity under different streamwise locations.

Figure 21 .
Figure 21.(a) The mean vorticity magnitude; (b) The mean streamwise velocity; (c) The mean horizontal velocity; (d) The mean vertical velocity.Comparison of the mean vorticity under different streamwise locations.

Figure 21 .
Figure 21.(a) The mean vorticity magnitude; (b) The mean streamwise velocity; (c) The mean horizontal velocity; (d) The mean vertical velocity.Comparison of the mean vorticity under different streamwise locations.J. Mar.Sci.Eng.2023, 11, x FOR PEER REVIEW 21 of 25

Figure 22 .
Figure 22.(a) The turbulence kinetic energy; (b) The normal stress in the streamwise component; (c) The normal stress in the horizontal component; (d) The normal stress in the vertical component.Comparison of the turbulence kinetic energy under different streamwise locations.

Figure 22 .
Figure 22.(a) The turbulence kinetic energy; (b) The normal stress in the streamwise component; (c) The normal stress in the horizontal component; (d) The normal stress in the vertical component.Comparison of the turbulence kinetic energy under different streamwise locations.

Table 1 .
Summary of the centers of the vortices on the upper hull.

Table 1 .
Summary of the centers of the vortices on the upper hull.

Table 2 .
Summary of the relative errors to Exp. 2 of the centers of the vortices.