Zhideng Zhou, Zhaobin Li, Guowei He, Xiaolei Yang
a The State Key Laboratory of Nonlinear Mechanics, Institute of Mechanics, Chinese Academy of Sciences, Beijing 100190, China
b School of Engineering Sciences, University of Chinese Academy of Sciences, Beijing 10 0 049, China
Keywords:Multi-fidelity model Large-eddy simulation Immersed boundary method Actuator disk model DARPA suboff
ABSTRACT Flow around a real-life underwater vehicle often happens at a high Reynolds number with flow structures at different scales from the boundary layer around a blade to that around the hull. This poses a great challenge for large-eddy simulation of an underwater vehicle aiming at resolving all relevant flow scales.In this work, we propose to model the hull with appendages using the immersed boundary method, and model the propeller using the actuator disk model without resolving the geometry of the blade. The proposed method is then applied to simulate the flow around Defense Advanced Research Projects Agency(DARPA) suboff. An overall acceptable agreement is obtained for the pressure and friction coefficients.Complex flow features are observed in the near wake of suboff. In the far wake, the core region is featured by a jet because of the actuator disk, surrounded by an annular region with velocity deficit due to the body of suboff.
Flows around the sterns of ships and underwater vehicles are complicated, with unsteady separations caused by the pressure gradients and the hull curvature, and are featured by flow structures at different scales from small protrusions on the surface, to propeller blade and the hull. The computational demand is extremely high for directly simulating all relevant scales. In this work, we propose to simulate the flow over the hull with appendages using the immersed boundary method, and model the propeller using the actuator disk model, respectively, and apply the proposed method to Defense Advanced Research Projects Agency(DARPA) suboff.
DARPA suboff is a standard model for investigating the hydrodynamics of underwater vehicles. Measurements on the flow over the suboff or suboff-like body have been carried out in the literature, for instance, the experiments by Huang et al. [1] and Jiménez et al. [2,3] , which provide useful dataset for validating computational models. Numerical simulations of suboff have also been carried out in the literature [4,5] . Kumar and Mahesh [6] employed wall-resolved large-eddy simulation (WRLES) to simulate the flow over the bare hull of suboff. Posa and Balaras [7] carried out WRLES of the suboff with appendages, and observed a bimodal behaviour for the turbulent stresses in the wake. Often the Reynolds number, which can be feasible for WRLES, is relatively low. Posa and Balaras [8] further examined the Reynolds number effects for suboff, found that the growth of boundary layer thickness is almost independent of the Reynolds number. To alleviate the grid resolution requirement near the surface, wall-modelled LES (WMLES) can be employed. Alin et al. [9] reviewed the WMLES, detached eddy simulation (DES) and Reynolds-averaged Navier-Stokes(RANS) for simulating the flow around the suboff with and without appendages. Liefvendahl and Fureby [10] reviewed the grid resolution requirements for WRLES and WMLES of ship hull hydrodynamics, and indicated that WMLES is the only feasible way for fullscale simulations. Grid generation is another challenge for simulating suboff with appendages and propeller. The immersed boundary(IB) method based on non-boundary-conforming grids can easily deal with complex geometries [11–13] . Simulation of suboff using the IB method can be found in the work by Shi et al. [14] .
The suboff is often simulated without a propeller. In the literature, the flow over a propeller has been simulated with different fidelity. Kumar and Mahesh [15] simulated the flow of a five-bladed propeller using WRLES. Liao et al. [16] simulated the flow around a propeller in crashback mode using WMLES with the curvilinear immersed boundary (CURVIB) method. To alleviate the grid need for geometry-resolved simulations, Liao et al. [17] proposed an actuator surface model for propeller, which employs a separate RANS simulation to compute the force coefficients on the surface,and apply the obtained force coefficients to the actuator surface for propeller in the LES of propeller wake. Simulation of a propeller using WRLES is already very expensive, for instance in the work by Balaras et al. [18] , 3.3 billion grid nodes were employed for simulating the INSEAN E1619 propeller under openwater condition. Simulating the suboff with a propeller is even more expensive. In the work by Posa and Balaras [19] , 3.5 billion nodes were employed for WRLES of suboff with appendages and propeller.
To simulate the flow over the suboff with a propeller in a computationally efficient way, we propose to model the propeller using the actuator disk model, which has been widely employed for wind turbine applications [20] , and simulate the flow over suboff body using the IB method with wall model, respectively, and test the proposed method by simulating the flow over the suboff with appendages and propeller under a typical condition. It is noted that the actuator disk model was also employed by Alin et al. [9] to model the effects of the propulsor, but with the bodyfitted method.
The proposed multi-fidelity method is composed of two parts,i.e., the IB method for suboff with appendages, and the actuator disk model for propeller. The proposed method is implemented in the virtual flow simulator (VFS-Wind) code [21] , which solves the incompressible Navier-Stokes equations in the curvilinear coordinate as follows:

The curvilinear immersed boundary method [23,24] is employed to simulate the hull of the underwater vehicle (e.g., the suboff body with appendages). For different variants of the CURVIB method developed for different applications, readers can refer to these papers [25–29] . In the CURVIB method, the background grid nodes are classified into fluid nodes and solid nodes. The fluid nodes with at least one neighbor in the solid are further identified as the IB nodes, at which the boundary conditions for the flow simulation are applied. For WRLES and DNS, the velocity at the IB nodes is linearly interpolated from the fluids nodes and the velocity at the boundary in the wall-normal direction. For WMLES, a simple linear interpolation is not suitable. The tangential velocity at the IB node is computed using the wall model, while the wall normal velocity is still computed using the linear interpolation. In the present study, the simplified thin boundary layer equation is solved in the wall model, which is written as

whereyis the coordinate in the wall-normal direction,uis the velocity tangential to the immersed boundary, andνtis the turbulent viscosity modelled using the mixing-length model with the van Driest damping function.
The actuator disk model is employed to take into account the effect of propeller on the flow. In the actuator disk model for wind turbine, the thrust is determined using the one-dimensional momentum theory [30] . In the present application, the force on the propeller F propeller , on the other hand, depends on the instantaneous hydrodynamic force F suboff on the body of suboff and the operational need of suboff (to accelerate, to decelerate, and etc.),which can be expressed as follows:

whereMsuboff and V suboff are the mass and the velocity of the suboff, respectively. The obtained force is then distributed on the actuator disk for propeller. In the present work, a relatively simple situation is considered, that (1) the acceleration of the suboff is zero; (2) the force distribution on the actuator disk is assumed being uniform; (3) only the hydrodynamic force in the axial direction is considered; (4) the tangential force exerted by a rotating propeller is not taken into account. At last, the obtained force on the actuator disk is spread to the surrounding background grid nodes using the smoothed discrete delta function [31] for advancing the flow to the next time step.
The setup for simulating the flow over the DARPA suboff is presented. The geometry of the DARPA suboff with appendages (AFF8)is shown in Fig. 1 . As seen, the bare hull of the suboff is composed of a streamlined forebody, a parallel middle body, and a stern with contraction in the radial direction; the appendages consist of the sail ( Fig. 1 b) and four stern fins ( Fig. 1 c). The propeller is modelled using the actuator disk model, and placed downstream of the fins as shown Fig. 1 c. Based on the size of the propeller in the literature [32] , the actuator disk of radius 0.258Dis placed at 8.225Dfrom the upstream edge of the suboff, whereDdenotes the diameter of the middle body. In the simulation, the Reynolds number based on the length of the suboffL= 8.6Dis equal toRe=U0L/ν= 1.2 ×107, whereU0 denotes the velocity of incoming flow,νdenotes the kinematic viscosity of fluid. Uniform inflow superimposed with velocity fluctuations generated by a von Kármán spectrum [33] (turbulence intensity: about 1% of the uniform flow)is applied at the inlet.
The suboff is placed with 0 angle of attack and 0 yaw angle,as shown in Fig. 2 . The sizes of the computational domain in the streamwise, normal and spanwise directions are [ ?2.6D, 23.2D] ×[ ?4.3D, 4.3D] ×[ ?4.3D, 4.3D], with grid numbers 841 ×361 ×361 ,respectively. Table 1 shows the details of grid spacing and number of grid nodes at different locations, and the maximum size of sail is [ 1.82D, 2.54D] ×[0.495D, 0.905D] ×[ ?0.0 6 6D, 0.0 6 6D],the maximum size of the fin in the positiveydirection is [ 7.49D,7.89D] ×[0.162D, 0.48D] ×[ ?0.041D, 0.041D].

Fig. 1. ( a ) Full geometry of the DARPA suboff with appendages; ( b ) sail on the parallel middle body; ( c ) fins and actuator disk for propeller.

Fig. 2. ( a ) Computational domain with the xOy slice of the Cartesian background grid and the suboff with appendages (1 in every 5 nodes shown); ( b ) the unstructured surface grid for the suboff body with sail; ( c ) the unstructured surface grid for the actuator disk and fins.

Fig. 3. Comparisons of ( a ) pressure coefficient C p and ( b ) skin-friction coefficient C f obtained from the present WMLES, the WRLES of Posa and Balaras [8] and experiment of Huang et al. [1] .

Table 1 Details of the grid spacing ( Δh ) and number of grid nodes ( N) at different locations,in which the coordinates and the grid spacing are normalized by D . In the row of grid spacing, the abbreviation “u” denotes the uniform grid, “r” and “l” denote the non-uniform grid defined using the Tanh function, with the smallest grid cell on the right and left, respectively.
Simulation results are presented and discussed. In Fig. 3 , we first compare the obtained pressure and friction coefficients with those in the literature. It is seen in Fig. 3 a that the computedCpin general agrees with those in the literature. One difference is observed on top of the sail starting fromx/L= 0.21 , where theCpcomputed in this work is lower than those in the literature. The other difference is located behind the actuator disk, which is located atx/L≈0.956 , theCpis higher than those without an actuator disk due to the additional pressure jump induced by the thrust on the actuator disk. The friction coefficientCfis compared with those in the literature in Fig. 3 b. Overall, the agreement is acceptable considering the relatively coarse grid employed in this work. The major differences are observed at locations upstream of the sail, where theCfpredicted by the present work is lower than those in the literature. One possible reason causing this difference is the numerical trip wire employed by Posa and Balaras [8] for triggering transition to turbulence, which is not employed in this work. The other difference is observed at the downstream locations of the sail, where theCfpredicted in this work is lower than those in the literature. This is related to how well the flow over the sail is predicted, which is affected by the grid resolution as well as the inflow, which, on the other hand, is influenced by the numerical trip wire placed close to the upstream end of suboff as in the work by Posa and Balaras [8] .

Fig. 4. Time-averaged streamwise velocity at ( a ) xOy slice ( z = 0 . 0 ) and ( b ) the slice rotating counterclockwise with the angle of 45 ?from xOy slice ( y + z = 0 ).

Fig. 5. Time-averaged pressure distribution at the two-dimensional slice of ( a ) z = 0 and ( b ) y + z = 0 .

Fig. 6. Time-averaged streamwise velocity U at the yOz slices with x/D = 1 ~10 .

Fig. 7. Contours of profiles of ( a ) time-averaged streamwise velocity, ( b ) time-averaged pressure and ( c ) variance of streamwise velocity fluctuations on the xOy plane located at z = 0 . The profiles are located at x/D = 8 . 7 , and 9 . 0 ~18 . 0 with increment of 1, respectively.
In order to show the impact of having an actuator disk model on the global flow features around the suboff with appendages,Figs. 4 and 5 depict the contours of time-averaged streamwise velocityUand pressurePat two slices along the centerline of the bare hull at different azimuthal angles, one withz= 0 passing through the sail and one fin (a), one withy+z= 0 without cutting any appendages (b). As seen (a) of both figures, the appendages significantly impact the flow around them, i.e., low speed wake region in the downstream, and high pressure region in the upstream.The flow composed of a jet region and a surrounding wake region is observed extending for a long distance. The width of the jet region is observed being higher in theydirection as shown in Fig. 4 a when compared with that shown in Fig. 4 b. As for the pressure, variations are large in the near wake region and decay to zero quickly as travelling downstream. Then we examine in Fig. 6 the flow field atyOzslices located at different streamwise locations. As seen, the flow atx/D= 1 ~7 is featured by boundary layer at most azimuthal locations, except in the downstream of the sail where its wake dominates. Atx/D= 8 , where the flow passes through the stern with contraction in the radial direction, the flow exhibits significant deceleration. Atx/D= 9,10 in the near wake of the suboff,the flow around the centerline of suboff is significantly accelerated due to the thrust on the disk. The last considered variation lies in the boundary layer at the parallel middle body. The boundary layer atz= 0,y> 0 is thicker than that aty+z= 0 because of the deceleration of near-wall velocity.
After showing the flowfield around the body and in the near wake. Here, we examine in Fig. 7 the flowfield in the far wake with profiles of flow statistics atx/D= 8.7 and 9.0,10.0,11.0,···,18.0 ,respectively. As seen, the time-averaged streamwise velocity and pressure, and the variance of streamwise velocity fluctuations all variate significantly in the near wake (untilx/D= 10 , i.e., 1.4Dfrom the downstream end of suboff) in the transverse directions.Atx/D= 8.7,9 locations, it is seen that the streamwise velocity profile is featured by velocity deficit in three banded regions, i.e.,one located along the centerline and the other two located at two ends, and two jet-like regions located in-between the three velocity deficit regions. As for the pressure, a peak is observed at these two locations with an abrupt decrease in its magnitude fromx/D= 8.7 to 9. The transverse variations of 〈u′u′〉 are observed being fairly complex at these two locations due to the complex geometry of suboff near its downstream end and the actuator disk.At further downstream locations starting fromx/D= 10.0 , the tranverse variations of different flow characteristics become relatively simple. For the streamwise velocity, the jet due to the actuator disk for propeller occupies most of the central area of the wake, surrounded by the region with velocity deficit due to the appendages of suboff. The transverse variations of the pressure, on the other hand, are fairly small at these far wake locations. As for the variance of the streamwise velocity fluctuations, a region with peaks is observed in the wake at these far wake locations.
In this work, we propose to simulate the flow around an underwater vehicle in a multi-fidelity way, that its hull with appendages is geometrically resolved using the immersed boundary method,while the propeller is parameterized using the actuator disk model,and then apply this strategy to simulate the flow over DARPA suboff with appendages and propeller.
To evaluate the simulation results, we compare the obtained pressure and friction coefficients with those in the literature. The pressure coefficient in general agrees with those in the literature.As for the friction coefficient, under-predictions are observed at locations upstream and downstream of the sail, which is probably due to the employed coarse grid, and the fact that a numerical trip wire located near upstream end of suboff for triggering turbulence transition is not employed in this work, which is often employed in the literature. The flowfield around the body and in the wake are then examined, with complicated flow structures observed near the downstream end of suboff. At further downstream locations, a wake composed of a core region of fast flowing fluid and a surrounding region of slow flowing fluid is observed extending for a long distance.
Further development is needed to apply the proposed method to simulate a full-scale underwater vehicle, e.g., wall models [34] ,which can accurately take into account the effects of the pressure gradient and other non-equilibrium effects, and parameterized models, which can accurately model the effect of small scale structures on the surface [35] .
Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.
Acknowledgment
This work was supported by the National Natural Science Foundation of China (NSFC) Basic Science Center Program for “Multiscale Problems in Nonlinear Mechanics” (No. 11988102), NSFC(No. 12002345), and China Postdoctoral Science Foundation (No.2020M680027).
Theoretical & Applied Mechanics Letters2022年1期