Quasi-Analytical Solution of Optimum and Maximum Depth of Transverse V-Groove for Drag Reduction at Different Reynolds Numbers

: Reducing the skin-friction drag of a vehicle is an important way to reduce carbon emissions. Previous studies have investigated the drag reduction mechanisms of transverse grooves. However, it is more practical to investigate which groove geometry is optimal at different inﬂow conditions for engineering. The purpose of this paper is to establish the physical model describing the relationship between the dimensionless depth ( H + = Hu τ / υ ) of the transverse groove, the dimensionless inﬂow velocity ( U + ∞ = U ∞ / u τ ), and the drag reduction rate ( η ) to quasi-analytically solve the optimal and maximum transverse groove depth according to the Reynolds numbers. Firstly, we use the LES with the dynamic subgrid model to investigate the drag reduction characteristics of transverse V-grooves with different depths (h = 0.05~0.9 mm) at different Reynolds numbers (1.09 × 10 4 ∼ 5.44 × 10 5 ) and ﬁnd that H + and U + ∞ affect the magnitude of slip velocity ( U + s ), thus driving the variation of the viscous drag reduction rate ( η ν ) and the increased rate of pressure drag ( η p ). Moreover, the relationship between U + s , η ν , and η p is established based on the slip theory and the law of pressure distribution. Finally, the quasi-analytical solutions for the optimal and maximum depths are solved by adjusting U + s to balance the cost ( η p ) and beneﬁt ( η ν ). This solution is in good agreement with the present numerical simulations and previous experimental results. over the grooved plate, respectively. It is apparent that the simulation results show good predictions of the location and structure of the vortex formed within the groove, which indicates that the CFD method can accurately simulate the flow details over the plate.


Introduction
In recent years, due to the increase in energy consumption and the strict requirements of fuel efficiency, the technologies of reducing carbon emissions associated with skinfriction drag reduction have drawn much attention [1]. For Lufthansa Cargo's Boeing 777F freighters, reducing the skin-friction drag by 1% means annual savings of around 3700 tons of kerosene and just under 11,700 tons of CO 2 emissions [2]. Compared with traditional drag reduction methods, bionic microstructures have a better potential for engineering applications because of their remarkable drag reduction properties and good applicability [3][4][5]. Previous studies have shown that there are two types of microstructure, one is the riblets imitating shark shin [6,7] and the other is the transverse grooves imitating dolphin skin [8][9][10][11], which are parallel and perpendicular to the flow direction, respectively. It has been reported that longitudinal riblets are capable of delivering a reduction of surface friction drag around 10% [12]. The drag reduction mechanism of longitudinal riblets is attributed to the damping of crossflow fluctuations or the uplift of turbulent streamwise vortices above the riblet valley [13][14][15].
The better drag reduction properties of transverse grooves have been proved by the latest research. Lee et al. [16] observed that the maximum measured drag reduction of 40% was achieved by a nanoporous transverse grooved plate. Liu et al. [17] used large eddy simulation (LES) technology to analyze entropy generation in the flow over a transverse grooved plate. The results showed that the total entropy generation in the near-wall region decreased by approximately 25%. Wang et al. [18] reported that an 18.76% net drag reduction was achieved by transverse grooves. These studies indicate that transverse grooves exhibit more significant potential for engineering applications than streamwise riblets.
At present, most studies focus on the drag reduction mechanisms of transverse grooves. Some studies suggested that the vortices formed within the grooves weaken the turbulence structure in the boundary layer near the wall [19,20]. The more popular perspective indicated that the vortices within the transverse grooves change the sliding friction into rolling friction at the solid-liquid interface, which is also named the "micro air-bearing phenomenon" [21][22][23][24][25]. Several studies focused on the Cassie-Baxter state to reduce the solidliquid contact area and create a superhydrophobic effect over the transverse grooves [26][27][28]. The studies of Seo et al. [29], Lang et al. [30], and Mariotti et al. [31] showed that the vortices formed in the transverse grooves increase the momentum in the boundary layer near the wall, thus effectively controlling the flow separation.
All the above studies are of great significance to qualitatively explain the drag reduction mechanism of the transverse groove. However, it is of more practical interest to understand which groove geometry is optimal at different inflow conditions for engineering. To the best of the authors' knowledge, compared with a large number of studies in the context of streamwise riblets, less is known regarding the parameter studies on transverse grooves. The drag reduction characteristics of transverse grooves are mainly determined by the shape, AR (ratio of groove width to depth), and depth. Cui et al. [22] conducted a numerical simulation on the pressure drop in microchannel flow over different transverse-grooved surfaces and found that the drag-reduction rate of V-shaped transverse grooves is better than that of rectangular transverse grooves. The experimental results of Liu et al. [32] demonstrated that the drag reduction performance was best when the AR is 2 at a Reynolds number of 50,000. The purpose of this paper is to establish the physical model describing the relationship between the dimensionless depth (H + = Hu τ /υ) of the transverse groove, the dimensionless inflow velocity (U + ∞ = U ∞ /u τ ), and the drag reduction rate (η), so as to quasi-analytically solve the optimal and maximum transverse groove depth according to the different Reynolds numbers. The optimum depth of the transverse groove corresponds to the maximum drag reduction rate, and the maximum depth corresponds to the drag reduction rate of zero, which is the limit for the allowable machining error for engineering applications.
This paper is organized as follows. First, the numerical methodology is formulated in Section 2. Secondly, the drag reduction characteristics of transverse V-grooves with different depths at different Reynolds numbers are investigated by LES in Section 3. Then, in Section 4, the theoretical model for the optimal and maximum depth of the transverse groove is established. Based on this model, the quasi-analytical solution of the groove depth has been solved and several grooved plates with different groove depths have been designed to verify the drag reduction characteristics in Section 5. Finally, the conclusions are presented in Section 6.

Solving Methods
The large eddy simulation (LES) method is used in the commercial software FLUENT 18.0 to obtain the induced drag reduction and flow characteristics [33,34]. A dynamic subgrid-scale (SGS) is chosen to model the unresolved small flow field scale motions [35,36]. The discretized continuity equation is solved using the Rhie and Chow method [37] to compute the mass flux at each face. The diffusion terms and the advection terms in the discretized momentum equations are solved using a second-order-accurate centraldifferencing discretization scheme and a second-order upwind scheme, respectively. Moreover, a second-order implicit time-stepping approach is used for temporal discretization. Thus, the space and time resolution of the numerical method is of second-order accuracy. For the pressure-velocity coupling, SIMPLE (Semi-Implicit method for pressure linked equations) is used to enforce the mass conservation and obtain the pressure fields. For the evaluation of gradients and derivatives, the least-squares cell-based gradient method is employed. The dimensionless physical timestep ∆tU/H ≈ 0.02 [38] is used, where U denotes the uniform velocity at the inlet and H represents the depth of grooves. The dimensionless time of statistical averaging is above TU/H ≈ 400 [34] to obtain the time-averaged results.

Computational Domain and Boundary Conditions
The 3D computational domain and the boundary conditions are shown schematically in Figure 1 and Table 1. The height and spanwise length of the computational domain are 25 mm (L + y = L + z = 103) [34], which covers at least the thickness of the boundary layer under all the flow conditions. The total length of the computational domain is 200 mm (L + x = 833). A smooth wall with a length of 160 mm is placed upstream of a grooved wall for the development of turbulence from laminar flow. Another smooth wall with a length of 20 mm is located downstream of the grooved wall to prevent the propagations of pressure perturbations at the outlet. The simulated grooved wall is about 20 mm long, consisting of different symmetric V-grooves profiles-whose AR is 2 (aspect ratios, S/H, where S represents the width of an individual groove) and depths are 0.05~0.9 mm (a baseline plate without grooves is also simulated in the same position). Reynolds numbers range from 1.09 × 10 4 to 5.44 × 10 5 , which is based on the length of the flat wall placed upstream of a grooved wall (160 mm), the freestream velocity, and the viscosity of the fluid. At all the solid walls, the no-slip condition is specified. At the inlet of the computational domain, an ideal gas flows with uniform velocity and a freestream turbulence intensity of 1%, and the pressure condition is set at the outlet. discretized momentum equations are solved using a second-order-accurate central-differencing discretization scheme and a second-order upwind scheme, respectively. Moreover, a second-order implicit time-stepping approach is used for temporal discretization. Thus, the space and time resolution of the numerical method is of second-order accuracy. For the pressure-velocity coupling, SIMPLE (Semi-Implicit method for pressure linked equations) is used to enforce the mass conservation and obtain the pressure fields. For the evaluation of gradients and derivatives, the least-squares cell-based gradient method is employed. The dimensionless physical timestep Δ / 0.02 [38] is used, where U denotes the uniform velocity at the inlet and H represents the depth of grooves. The dimensionless time of statistical averaging is above / 400 [34] to obtain the time-averaged results.

Computational Domain and Boundary Conditions
The 3D computational domain and the boundary conditions are shown schematically in Figure 1 and Table 1. The height and spanwise length of the computational domain are 25 mm ( = = 103) [34], which covers at least the thickness of the boundary layer under all the flow conditions. The total length of the computational domain is 200 mm ( = 833). A smooth wall with a length of 160 mm is placed upstream of a grooved wall for the development of turbulence from laminar flow. Another smooth wall with a length of 20 mm is located downstream of the grooved wall to prevent the propagations of pressure perturbations at the outlet. The simulated grooved wall is about 20 mm long, consisting of different symmetric V-grooves profiles-whose AR is 2 (aspect ratios, S/H, where S represents the width of an individual groove) and depths are 0.05~0.9 mm (a baseline plate without grooves is also simulated in the same position). Reynolds numbers range from 1.09 × 10 to5.44 × 10 , which is based on the length of the flat wall placed upstream of a grooved wall (160 mm), the freestream velocity, and the viscosity of the fluid. At all the solid walls, the no-slip condition is specified. At the inlet of the computational domain, an ideal gas flows with uniform velocity and a freestream turbulence intensity of 1%, and the pressure condition is set at the outlet.    In order to verify that the normal height (L + y ) and spanwise length (L + z ) of the computational domain meet the requirements of LES, two normal heights (L + y = 103 and 123.6) and two spanwise lengths (L + z = 103 and 123.6) are chosen to investigate the effect on outcomes. Table 2 shows the simulation results of the total drag of the grooved plate. The drag hardly changes when the normal height and the spanwise length are greater than 103, which means that the results conducted with L + x = 103 and L + z = 103 in the present are adequate for the requirements of LES.

Grid Independence Study
The structured mesh is generated by commercial software, ICEM, as shown in Figure 2. The grid resolution and the number of grid nodes are shown in Table 1. The grids are clustered near the wall surface and the normal distance from the wall surface to the nearest grid points Y + is 0.04. The maximum normal grid resolution ∆y + max is less than 10. The streamwise grid resolution ∆x + is 0.07 within the grooves, and ∆x + max is less than 10 in other streamwise positions. The spanwise grid resolution ∆z + is 4.1 [34]. In order to verify that the normal height ( ) and spanwise length ( ) of the computational domain meet the requirements of LES, two normal heights ( = 103 and 123.6) and two spanwise lengths ( = 103 and 123.6) are chosen to investigate the effect on outcomes. Table 2 shows the simulation results of the total drag of the grooved plate. The drag hardly changes when the normal height and the spanwise length are greater than 103, which means that the results conducted with = 103 and = 103 in the present are adequate for the requirements of LES.

Grid Independence Study
The structured mesh is generated by commercial software, ICEM, as shown in Figure  2. The grid resolution and the number of grid nodes are shown in Table 1. The grids are clustered near the wall surface and the normal distance from the wall surface to the nearest grid points Y is 0.04. The maximum normal grid resolution ∆ is less than 10. The streamwise grid resolution ∆ is 0.07 within the grooves, and ∆ is less than 10 in other streamwise positions. The spanwise grid resolution ∆ is 4.1 [34]. In order to verify that the grid resolution meets the requirements of the large eddy simulation, two streamwise grid resolutions within the groove ∆ (0.07 and 0.32) and two spanwise grid resolutions ∆ (4.1 and 10) are chosen to investigate the effect on the outcomes. Table 3 shows the simulation results for the drag of the grooved plate and the streamline inside the groove. The resistances of a grooved plate hardly change when ∆ = 0.07 (groove) and ∆ = 4.1, which are selected for deriving all the other results.  In order to verify that the grid resolution meets the requirements of the large eddy simulation, two streamwise grid resolutions within the groove ∆x + (0.07 and 0.32) and two spanwise grid resolutions ∆z + (4.1 and 10) are chosen to investigate the effect on the outcomes. Table 3 shows the simulation results for the drag of the grooved plate and the streamline inside the groove. The resistances of a grooved plate hardly change when ∆x + = 0.07 (groove) and ∆z + = 4.1, which are selected for deriving all the other results. In order to verify that the normal height ( ) and spanwise length ( ) of the computational domain meet the requirements of LES, two normal heights ( = 103 and 123.6) and two spanwise lengths ( = 103 and 123.6) are chosen to investigate the effect on outcomes. Table 2 shows the simulation results of the total drag of the grooved plate. The drag hardly changes when the normal height and the spanwise length are greater than 103, which means that the results conducted with = 103 and = 103 in the present are adequate for the requirements of LES.

Grid Independence Study
The structured mesh is generated by commercial software, ICEM, as shown in Figure  2. The grid resolution and the number of grid nodes are shown in Table 1. The grids are clustered near the wall surface and the normal distance from the wall surface to the nearest grid points Y is 0.04. The maximum normal grid resolution ∆ is less than 10. The streamwise grid resolution ∆ is 0.07 within the grooves, and ∆ is less than 10 in other streamwise positions. The spanwise grid resolution ∆ is 4.1 [34]. In order to verify that the grid resolution meets the requirements of the large eddy simulation, two streamwise grid resolutions within the groove ∆ (0.07 and 0.32) and two spanwise grid resolutions ∆ (4.1 and 10) are chosen to investigate the effect on the outcomes. Table 3 shows the simulation results for the drag of the grooved plate and the streamline inside the groove. The resistances of a grooved plate hardly change when ∆ = 0.07 (groove) and ∆ = 4.1, which are selected for deriving all the other results. In order to verify that the normal height ( ) and spanwise length ( ) of the computational domain meet the requirements of LES, two normal heights ( = 103 and 123.6) and two spanwise lengths ( = 103 and 123.6) are chosen to investigate the effect on outcomes. Table 2 shows the simulation results of the total drag of the grooved plate. The drag hardly changes when the normal height and the spanwise length are greater than 103, which means that the results conducted with = 103 and = 103 in the present are adequate for the requirements of LES.

Grid Independence Study
The structured mesh is generated by commercial software, ICEM, as shown in Figure  2. The grid resolution and the number of grid nodes are shown in Table 1. The grids are clustered near the wall surface and the normal distance from the wall surface to the nearest grid points Y is 0.04. The maximum normal grid resolution ∆ is less than 10. The streamwise grid resolution ∆ is 0.07 within the grooves, and ∆ is less than 10 in other streamwise positions. The spanwise grid resolution ∆ is 4.1 [34]. In order to verify that the grid resolution meets the requirements of the large eddy simulation, two streamwise grid resolutions within the groove ∆ (0.07 and 0.32) and two spanwise grid resolutions ∆ (4.1 and 10) are chosen to investigate the effect on the outcomes. Table 3 shows the simulation results for the drag of the grooved plate and the streamline inside the groove. The resistances of a grooved plate hardly change when ∆ = 0.07 (groove) and ∆ = 4.1, which are selected for deriving all the other results.  Figure 3 shows the relative error of the total resistance of a grooved plate (with symmetric V-grooves whose AR is 2) compared with that of a smooth plate (based on the total drag in the case of 2,403,465 grids) at five different grid-refinement levels. The resistances of a grooved plate and the baseline plate no longer change when the number of grid cells is greater than 1,603,491, and the relative difference is 0.0465%, which is used as the grid resolution in deriving all the other results.
Symmetry 2022, 14, 342 5 of 25 Figure 3 shows the relative error of the total resistance of a grooved plate (with symmetric V-grooves whose AR is 2) compared with that of a smooth plate (based on the total drag in the case of 2,403,465 grids) at five different grid-refinement levels. The resistances of a grooved plate and the baseline plate no longer change when the number of grid cells is greater than 1,603,491, and the relative difference is 0.0465%, which is used as the grid resolution in deriving all the other results.

Experimental Validation
In order to check the accuracy of the obtained numerical results, the results of the numerical simulation are compared with the relevant experimental data reported previously. The numerical results of velocity profiles over the grooved plates with a groove depth of 1.62 mm and groove width of 3.57 mm are compared with the experimental results obtained by Ahmadi-Baloutaki et al. [39] in Figure 4a,b (the Reynolds number based on the length of the grooved plate is 1.85 × 10 , and the turbulence intensity is 0.4% and 4.4%, respectively). The simulated velocity profiles agree well with the previously published experimental results, suggesting that the present numerical approach has sufficient accuracy for predicting the drag reduction of transverse grooves (the relative error is found to be less than 3%). Moreover, Figure 4c,d show the experimental results [40] and numerical results of the velocity vector over the grooved plate, respectively. It is apparent that the simulation results show good predictions of the location and structure of the vortex formed within the groove, which indicates that the CFD method can accurately simulate the flow details over the grooved plate.

Experimental Validation
In order to check the accuracy of the obtained numerical results, the results of the numerical simulation are compared with the relevant experimental data reported previously. The numerical results of velocity profiles over the grooved plates with a groove depth of 1.62 mm and groove width of 3.57 mm are compared with the experimental results obtained by Ahmadi-Baloutaki et al. [39] in Figure 4a,b (the Reynolds number based on the length of the grooved plate is 1.85 × 10 5 , and the turbulence intensity is 0.4% and 4.4%, respectively). The simulated velocity profiles agree well with the previously published experimental results, suggesting that the present numerical approach has sufficient accuracy for predicting the drag reduction of transverse grooves (the relative error is found to be less than 3%). Moreover, Figure 4c,d show the experimental results [40] and numerical results of the velocity vector over the grooved plate, respectively. It is apparent that the simulation results show good predictions of the location and structure of the vortex formed within the groove, which indicates that the CFD method can accurately simulate the flow details over the grooved plate.
Symmetry 2022, 14, 342 5 of 25 Figure 3 shows the relative error of the total resistance of a grooved plate (with symmetric V-grooves whose AR is 2) compared with that of a smooth plate (based on the total drag in the case of 2,403,465 grids) at five different grid-refinement levels. The resistances of a grooved plate and the baseline plate no longer change when the number of grid cells is greater than 1,603,491, and the relative difference is 0.0465%, which is used as the grid resolution in deriving all the other results.

Experimental Validation
In order to check the accuracy of the obtained numerical results, the results of the numerical simulation are compared with the relevant experimental data reported previously. The numerical results of velocity profiles over the grooved plates with a groove depth of 1.62 mm and groove width of 3.57 mm are compared with the experimental results obtained by Ahmadi-Baloutaki et al. [39] in Figure 4a,b (the Reynolds number based on the length of the grooved plate is 1.85 × 10 , and the turbulence intensity is 0.4% and 4.4%, respectively). The simulated velocity profiles agree well with the previously published experimental results, suggesting that the present numerical approach has sufficient accuracy for predicting the drag reduction of transverse grooves (the relative error is found to be less than 3%). Moreover, Figure 4c,d show the experimental results [40] and numerical results of the velocity vector over the grooved plate, respectively. It is apparent that the simulation results show good predictions of the location and structure of the vortex formed within the groove, which indicates that the CFD method can accurately simulate the flow details over the grooved plate.

Characteristics of Drag Reduction Induced by Transverse V-Grooves with Different Depths at Different Reynolds Numbers
Before establishing the prediction model of the optimal and maximum depth, it is necessary to analyze the key physical factors affecting the drag reduction characteristics. In this section, the amounts of drag reduction induced by grooves with different depths at different Reynolds numbers are computed and the mechanisms of drag reduction are analyzed. Figure 5 illustrates the variation in drag-reduction rate with depth at different Reynolds numbers. The drag-reduction rate is defined as,

=
( in which and represent the resistance of the grooved plate and the baseline plate, respectively. The results show that the drag-reduction rate induced by the grooved plate first increases and then decreases with the increase of depth. At each Reynolds number, there is an optimal groove depth for maximum drag reduction and a maximum groove depth corresponding to the zero-drag reduction rate, which is qualitatively consistent with the results of the previous study in the context of streamwise riblets [12]. Interestingly, the optimal and maximum depth decrease with the increase in Reynolds number, as shown in Table 4. It is worth noting that because the groove depth in the simulation cases is discrete, the optimal groove depth chosen from Table 4 actually corresponds to the 'near-maximum' drag reduction rate, and the maximum depth actually corresponds to the 'near-zero' drag reduction rate.

Characteristics of Drag Reduction Induced by Transverse V-Grooves with Different Depths at Different Reynolds Numbers
Before establishing the prediction model of the optimal and maximum depth, it is necessary to analyze the key physical factors affecting the drag reduction characteristics. In this section, the amounts of drag reduction induced by grooves with different depths at different Reynolds numbers are computed and the mechanisms of drag reduction are analyzed. Figure 5 illustrates the variation in drag-reduction rate with depth at different Reynolds numbers. The drag-reduction rate is defined as, in which F G and F R represent the resistance of the grooved plate and the baseline plate, respectively. The results show that the drag-reduction rate induced by the grooved plate first increases and then decreases with the increase of depth. At each Reynolds number, there is an optimal groove depth for maximum drag reduction and a maximum groove depth corresponding to the zero-drag reduction rate, which is qualitatively consistent with the results of the previous study in the context of streamwise riblets [12]. Interestingly, the optimal and maximum depth decrease with the increase in Reynolds number, as shown in Table 4. It is worth noting that because the groove depth in the simulation cases is discrete, the optimal groove depth chosen from Table 4 actually corresponds to the 'near-maximum' drag reduction rate, and the maximum depth actually corresponds to the 'near-zero' drag reduction rate.

Characteristics of Drag Reduction Induced by Transverse V-Grooves with Different Depths at Different Reynolds Numbers
Before establishing the prediction model of the optimal and maximum depth, it is necessary to analyze the key physical factors affecting the drag reduction characteristics. In this section, the amounts of drag reduction induced by grooves with different depths at different Reynolds numbers are computed and the mechanisms of drag reduction are analyzed. Figure 5 illustrates the variation in drag-reduction rate with depth at different Reynolds numbers. The drag-reduction rate is defined as,

=
( in which and represent the resistance of the grooved plate and the baseline plate, respectively. The results show that the drag-reduction rate induced by the grooved plate first increases and then decreases with the increase of depth. At each Reynolds number, there is an optimal groove depth for maximum drag reduction and a maximum groove depth corresponding to the zero-drag reduction rate, which is qualitatively consistent with the results of the previous study in the context of streamwise riblets [12]. Interestingly, the optimal and maximum depth decrease with the increase in Reynolds number, as shown in Table 4. It is worth noting that because the groove depth in the simulation cases is discrete, the optimal groove depth chosen from Table 4 actually corresponds to the 'near-maximum' drag reduction rate, and the maximum depth actually corresponds to the 'near-zero' drag reduction rate.    (2) [41]) to analyze the physical mechanism driving the variation in total drag.
here, u τ represents the local friction velocity.   (2) [41]) to analyze the physical mechanism driving the variation in total drag.
here, represents the local friction velocity.   The streamline patterns and velocity that contour over the grooved plate are shown in Figure 6. Some vortices formed in the grooves can be distinctly observed as perpendicular to the flow direction. These boundary vortices act as "air bearings", which separate the boundary layer from the solid wall, resulting in fluid sliding over the grooved plate. In order to measure the sliding degree, the dimensionless velocity profiles ( = / , and = / ) at the centerline of the grooved plate and the baseline plate at the corresponding position are compared (Figure 6, right). The results show that the velocity gradient over the grooved plate is less than that on the baseline plate, and an induced slip velocity ( ) on the horizontal line can be selected as the quantitative parameter to describe the slip phenomenon. Moreover, the comparison results of the different numerical cases indicate that the slip velocity is affected by the groove depth (h) and inflow velocity ( ).
The wall shear distribution and pressure contour shown in Figure 7 further illustrate that the magnitude of slip velocities drives the variation in the total viscous drag and pressure drag. On the one hand, these slip velocities reduce the velocity gradient over the groove plate, thus reducing the total viscous resistance. This point is further proved by the shear stress distribution diagram, which shows that the shear stress of the grooved The streamline patterns and velocity that contour over the grooved plate are shown in Figure 6. Some vortices formed in the grooves can be distinctly observed as perpendicular to the flow direction. These boundary vortices act as "air bearings", which separate the boundary layer from the solid wall, resulting in fluid sliding over the grooved plate. In order to measure the sliding degree, the dimensionless velocity profiles (Y + = yu τ /ν, and U + = U/u τ ) at the centerline of the grooved plate and the baseline plate at the corresponding position are compared (Figure 6, right). The results show that the velocity gradient over the grooved plate is less than that on the baseline plate, and an induced slip velocity (U + s ) on the horizontal line can be selected as the quantitative parameter to describe the slip phenomenon. Moreover, the comparison results of the different numerical cases indicate that the slip velocity is affected by the groove depth (h) and inflow velocity (U + ∞ ). The wall shear distribution and pressure contour shown in Figure 7 further illustrate that the magnitude of slip velocities drives the variation in the total viscous drag and pressure drag. On the one hand, these slip velocities reduce the velocity gradient over the groove plate, thus reducing the total viscous resistance. This point is further proved by the shear stress distribution diagram, which shows that the shear stress of the grooved wall is significantly less than that of the baseline plate. By comparing the numerical cases of grooved plates with the same depth at different inflow velocities, it can be seen that the larger the slip velocity, the smaller the velocity gradient. As well, the corresponding shear stress decreases the most, so the total viscous resistance decreases the most in turn. On the other hand, the "slip fluids" induced by the vortices separate on the leeward side of the groove and stagnation occurs on the windward side, resulting in additional pressure drag compared to the baseline plate-which increases the total resistance of the grooved plate. The larger the slip velocity, the greater the stagnation pressure on the windward side, resulting in greater additional pressure drag. In summary, the grooved plate reduces the viscous drag (benefits) and increases the pressure drag (costs), and the optimal drag reduction is the result of balancing the benefits and costs.
From the above analysis, the total drag of the grooved plate consists of the viscous drag (F GV ) and pressure drag (F GP ), which are expressed by Equations (3) and (4), respectively. F GV and F GP are determined by calculating the corresponding local stress-namely, shear (τ) and pressure (p) at the wall-and integrating the projected stress in the drag direction (e x , that is the unit vectors in the x direction shown in Figure 1 along the wetted wall (l s ).
here, p ∞ represents the ambient pressure, l is the unit area along the groove wall, and n denotes the normal vector to the wall. Therefore, Equation (1) is transformed into Equation (5): here, η ν = (F GV − F R )/F R denotes the reduction rate for viscous drag, and "η p = F GP /F R " denotes the increased rate of pressure drag. Figure 8 shows the variation of η ν and η p with U s , u t , and U + s (in order to compare the model results with the experimental results in the following sections, four different grooves with h of 0.05, 0.1, 0.15, and 0.3 mm are used for analysis). It can be seen that the absolute values of η ν and η p increase with the increase of U + s (U + s = U s /u t , where u t = τ ω /ρ is only a unit commonly used for dimensionless values. Therefore, there is no physical relationship between u t and the drag reduction rate), which further indicates that the slip velocity drives the variation in the total viscous drag and pressure drag. In conclusion, the essence of balancing the benefit (η ν ) and cost (η p ) to obtain the optimal drag reduction is how to obtain the optimal slip velocity by matching the depth of the groove and the inflow velocity.
Based on the above analysis, we propose a model to match the relationship between inflow velocity (U + ∞ ) and groove depth (h) by determining the optimal slip velocity (U + s ). The establishment process of this model is shown in Figure 9.
Step 3. Construct the physical relation between slip velocity (U + s ) and pressure dragincrease rate (η p ).
the absolute values of and increase with the increase of ( = / , where = ⁄ is only a unit commonly used for dimensionless values. Therefore, there is no physical relationship between and the drag reduction rate), which further indicates that the slip velocity d Based on the above analysis, we propose a model to match the relationship between inflow velocity ( ) and groove depth (h) by determining the optimal slip velocity ( ). The establishment process of this model is shown in Figure 9.
Step 2. Construct the physical relation between slip velocity ( ) and viscous drag-reduction rate ( ).
Step 3. Construct the physical relation between slip velocity ( ) and pressure drag-increase rate ( ).
Step 4. Balance and to design the groove depth ( ) in different inflow velocities ( ). Based on the above analysis, we propose a model to match the relationship between inflow velocity ( ) and groove depth (h) by determining the optimal slip velocity ( ). The establishment process of this model is shown in Figure 9.
Step 2. Construct the physical relation between slip velocity ( ) and viscous drag-reduction rate ( ).
Step 3. Construct the physical relation between slip velocity ( ) and pressure drag-increase rate ( ).
Step 4. Balance and to design the groove depth ( ) in different inflow velocities ( ).

Construction the Theoretical Model Describing the Relationship between the Dimensionless Depth (H + ), the Dimensionless Inflow Velocity (U + ∞ ), and Drag Reduction Rate (η)
4.1. The Relationship between Slip Velocity (U + s ), Depth (H + ), and Inflow Velocity (U + ∞ ) The streamline patterns shown in Figure 7 reveal that the boundary vortices are not full inside the grooves. Therefore, we assume that the distances from the vortex center to the groove bottom and the slip surface are Kh and L s (slip length, L s = (1 − K)h), respectively-as shown in Figure 10. Since the slip surface shown in Figure 10 (i.e., the position of the baseline plate) is at the viscous bottom layer, the total viscous stress can be expressed as

The Relationship between Slip Velocity ( ), Depth ( ), and Inflow Velocity ( )
The streamline patterns shown in Figure 7 reveal that the boundary vortices are not full inside the grooves. Therefore, we assume that the distances from the vortex center to the groove bottom and the slip surface are ℎ and (slip length, = 1 ℎ), respectively-as shown in Figure 10. Since the slip surface shown in Figure 10 (i.e., the position of the baseline plate) is at the viscous bottom layer, the total viscous stress can be expressed as Over the grooved plate, the velocity at the vortex core and the slip surface increases from 0 to , and both and are small quantities, so the velocity gradient of U in the Y direction at the midpoint of the slip surface can be approximated as = (7) in which = , then Equation (7) can be transformed into Equation (8). = therefore, substituting " = 1 ℎ" into Equation (8) yields the relation between the dimensionless slip velocity and dimensionless depth, which is described as here, is equal to / , is equal to ℎ/ , and K represents the fullness of the vortex within the groove, which varies with different inflow velocities and groove depths, as shown in the streamline patterns of Figure 11. Over the grooved plate, the velocity at the vortex core and the slip surface increases from 0 to U s , and both U s and L s are small quantities, so the velocity gradient of U in the Y direction at the midpoint of the slip surface can be approximated as ∂u ∂y in which τ ω = ρ(u τ ) 2 , then Equation (7) can be transformed into Equation (8).
therefore, substituting "L s = (1 − K)h" into Equation (8) yields the relation between the dimensionless slip velocity and dimensionless depth, which is described as here, U + s is equal to U s /u τ , H + is equal to ρu τ h/µ, and K represents the fullness of the vortex within the groove, which varies with different inflow velocities and groove depths, as shown in the streamline patterns of Figure 11.  To quantitatively estimate the variation of K, we use the Boltzmann function to fit the data points (see Figure 12), which can be expressed as Equation (10). The resulting fits are accurate to better than 2% at a 95% confidence interval.
where , , , and are the control parameters of the Boltzmann function, which can be described by Equation (11) (h is in millimeters).
therefore, substituting Equations (10) and (11) into Equation (9) yields Equation (12), which can predict the slip velocity over the groove plate with depth at the inflow velocity . Figure 13 compares the CFD data with the results of Equation (12). The results show that the prediction value of Equation (12) are qualitatively in agreement with the numerical simulations. In detail, when ℎ 0.1, the prediction accuracy decreases with the increase of , so Equation (12) is more adequate for the 14.88 < < 18.53. However, Figure 11. The position of boundary vortex varies with dimensionless velocity.
To quantitatively estimate the variation of K, we use the Boltzmann function to fit the data points (see Figure 12), which can be expressed as Equation (10). The resulting fits are accurate to better than 2% at a 95% confidence interval.
where A 1 , A 2 , u 0 , and du are the control parameters of the Boltzmann function, which can be described by Equation (11) (h is in millimeters).
12162h 2 U 0 = 18.5287 − 9.83878h + 6.70225h 2 du = 0.33396 + 4.79502h − 10.6357h 2 (11) therefore, substituting Equations (10) and (11) into Equation (9) yields Equation (12), which can predict the slip velocity over the groove plate with depth H + at the inflow velocity U + ∞ .  To quantitatively estimate the variation of K, we use the Boltzmann functi data points (see Figure 12), which can be expressed as Equation (10). The resul accurate to better than 2% at a 95% confidence interval.    Figure 13 compares the CFD data with the results of Equation (12). The results show that the prediction value of Equation (12) are qualitatively in agreement with the numerical simulations. In detail, when h ≤ 0.1, the prediction accuracy decreases with the increase of U + ∞ , so Equation (12) is more adequate for the 14.88 < U + ∞ < 18.53. However, when h > 0.1, Equation (12) has a greater prediction accuracy when the dimensionless inflow velocity is between 14.9 and 22.

The Relationship between Slip Velocity ( ) and Viscous Drag-Reduction Rate ( )
Fukagata and Kasagi [42] proposed a theoretical prediction model for the drag reduction rate achieved by superhydrophobic surfaces in a turbulent channel flow. By comparing the bulk mean velocity of the baseline plate and the hydrophobic surface, they successfully established the relationship between the streamwise and spanwise slip length and the drag reduction rate. Different from their work, our work aims to establish the relationship between slip velocity, , and viscous drag reduction rate, , over the grooved plate.
According to Equation (5), the viscous drag reduction rate is expressed by ⁄ , where ∝ and ∝ , so can be expressed as Equation (13).
here, and represent the local friction velocities of the grooved plate and baseline plate, respectively. Therefore, can be estimated if the relationship between and ⁄ is established.
The velocity profile over the grooved plate, as shown in Figure 14, can be regarded as the superposition of slip velocity, , and the no-slip flow velocity profile ( ), which can be described as Equation (14). Figure 15 < 1000), which was also demonstrated by the work of Min and Kim [44]. Therefore, and at = δ

The Relationship between Slip
Velocity (U + s ) and Viscous Drag-Reduction Rate (η ν ) Fukagata and Kasagi [42] proposed a theoretical prediction model for the drag reduction rate achieved by superhydrophobic surfaces in a turbulent channel flow. By comparing the bulk mean velocity of the baseline plate and the hydrophobic surface, they successfully established the relationship between the streamwise and spanwise slip length and the drag reduction rate. Different from their work, our work aims to establish the relationship between slip velocity, U + s , and viscous drag reduction rate, η ν , over the grooved plate. According to Equation (5), the viscous drag reduction rate η ν is expressed by (F GV − F R )/F R , where F GV ∝ u τg and F R ∝ u τb , so η ν can be expressed as Equation (13).
here, u τg and u τb represent the local friction velocities of the grooved plate and baseline plate, respectively. Therefore, η ν can be estimated if the relationship between U + s and u τg /u τb is established.
The velocity profile over the grooved plate, as shown in Figure 14, can be regarded as the superposition of slip velocity, U s , and the no-slip flow velocity profile (U g−s ), which can be described as Equation (14). Figure 15 shows the dimensionless velocity profiles of U g and U g−s (U + ∞ = 22, h = 0.1). The results show that if the local friction velocities of the grooved plate and baseline plate (u τg and u τb ) are used as the dimensionless units of U s and U g−s , respectively, then both U g /u τg and U g−s /u τb satisfy the logarithmic law (U + = 1 κ lnY + + B, κ = 0.41, B = 5.0) [43] in the logarithmic region (50 < Y + < 1000), which was also demonstrated by the work of Min and Kim [44]. Therefore, U g and U g−s at Y = K δ δ can be described as Equations (15) and (16), respectively, where δ stands for the boundary layer thickness and K δ is an adjustable constant (K δ = 0.1∼0.2 in the logarithmic region [41]). For the predicted equation, u τg over the grooved plate is unknown before numerical calculation or experimental measurement, so all u τg in Equation (15) is converted to can be described as Equations (15) and (16), respectively, where δ stands for the boundary layer thickness and is an adjustable constant ( = 0.1~0.2 in the logarithmic region [41]). For the predicted equation, over the grooved plate is unknown before numerical calculation or experimental measurement, so all in Equation (15) is converted to .  Similarly, in Equation (14) can be expressed as Equation (17) with .
With different in Figure 13 and the corresponding ∈ 14.9, 22.0 as an input.
Implicit Equation (18) is solved by the dichotomy to obtain . In Figure 16, we compare the prediction of Equation (18) with CFD data in the present, The CFD data reported by can be described as Equations (15) and (16), respectively, where δ stands for the boundary layer thickness and is an adjustable constant ( = 0.1~0.2 in the logarithmic region [41]). For the predicted equation, over the grooved plate is unknown before numerical calculation or experimental measurement, so all in Equation (15) is converted to .  Similarly, in Equation (14) can be expressed as Equation (17) with .
With different in Figure 13 and the corresponding ∈ 14.9, 22.0 as an input.
Implicit Equation (18) is solved by the dichotomy to obtain . In Figure 16, we compare the prediction of Equation (18) with CFD data in the present, The CFD data reported by Similarly, U s in Equation (14) can be expressed as Equation (17) with u τb .
Therefore, substituting Equations (13) and (15)- (17) into Equation (14) yields the relationship between η ν and U + s . With different U + s in Figure 13 and the corresponding U + ∞ ∈ [14.9, 22.0] as an input. Implicit Equation (18) is solved by the dichotomy to obtain η ν . In Figure 16, we compare the prediction of Equation (18) with CFD data in the present, The CFD data reported by Wu et al. (2019) [5], and the experimental data reported by Liu et al. (2020) [32]. The reason why the relationship between η ν and U + ∞ is used, η ν U + s , is shown in Figure 16 and is mainly to facilitate the comparison between the results predicted by Equation (18) and the data obtained by previous studies. It is worth noting that the Reynolds number in all the U + ∞ according to Equation (2), and the groove depths in our simulation are consistent with those in the references (h is 0.05, 0.1, and 0.3, respectively). 14.9 < U + ∞ < 20.35, the predicted drag reduction rates are in good agreement with those from the present numerical simulations and the previous experimental data U + ∞ > 20.35, there is a significant error between the prediction of the model and the numerical results of Wu et al. [5], which indicates that this model is not applicable at a high Reynolds number.  [5], and the experimental data reported by Liu et al. (2020) [32]. The reason why the relationship between and is used, ￼ ￼, is shown in Figure 16 and is mainly to facilitate the comparison between the results predicted by Equation (18) and the data obtained by previous studies. It is worth noting that the Reynolds number in all the ￼ according to Equation (2), and the groove depths in our simulation are consistent with those in the references (h is 0.05, 0.1, and 0.3, respectively). 14.9 < < 20.35￼, the predicted drag reduction rates are in good agreement with those from the present numerical simulations and the previous experimental data￼ 20.35￼, there is a significant error between the prediction of the model and the numerical results of Wu et al. [5], which indicates that this model is not applicable at a high Reynolds number.

The Relationship between Slip Velocity ( ) and Pressure Drag-Increase Rate ( )
Predicting the viscous drag reduction rate ( ) of the grooved plate according to the corresponding slip velocity ( ) has been conducted in Section 4.2. In this section, the relationship between and the pressure drag-increase rate ( ) will be established by analyzing the effect of the slip velocity on the pressure distribution of the groove wall.
According to Equation (5), the pressure drag-increase rate is expressed by Equation (19).
= (19) here, the drag of the baseline plate can be estimated as 20) where N denotes the number of transverse grooves in the corresponding length groove plate, and represents the dimensionless width of a groove. Figure 17 shows the gauge pressure distribution on the groove wall (h = 0.05) at different inflow velocities. The results show that the high-pressure region is formed with the stagnation of slip velocity induced by the boundary vortex on the groove windward side and the corresponding low-pressure region is formed on the groove leeward side due to

The Relationship between Slip Velocity (U +
s ) and Pressure Drag-Increase Rate (η p ) Predicting the viscous drag reduction rate (η ν ) of the grooved plate according to the corresponding slip velocity (U + s ) has been conducted in Section 4.2. In this section, the relationship between U + s and the pressure drag-increase rate (η p ) will be established by analyzing the effect of the slip velocity on the pressure distribution of the groove wall.
According to Equation (5), the pressure drag-increase rate η p is expressed by Equation (19).
here, the drag of the baseline plate can be estimated as where N denotes the number of transverse grooves in the corresponding length groove plate, and S + represents the dimensionless width of a groove. Figure 17 shows the gauge pressure distribution on the groove wall (h = 0.05) at different inflow velocities. The results show that the high-pressure region is formed with the stagnation of slip velocity induced by the boundary vortex on the groove windward side and the corresponding low-pressure region is formed on the groove leeward side due to local flow separation. Therefore, the additional pressure drag (F GP ) of the grooved plate is equal to the integral of the difference between the high and low pressure in each groove. For a groove, taking the axis with zero-gauge pressure as the abscissa axis and the centerline of the groove as the ordinate axis, the pressure distribution is a function of x, as shown in Figure 18. Therefore, the total pressure drag of the grooved plate is expressed as here θ is 45 degrees (AR = 2).
Symmetry 2022, 14, 342 16 of 25 local flow separation. Therefore, the additional pressure drag ( ) of the grooved plate is equal to the integral of the difference between the high and low pressure in each groove. For a groove, taking the axis with zero-gauge pressure as the abscissa axis and the centerline of the groove as the ordinate axis, the pressure distribution is a function of x, as shown in Figure 18. Therefore, the total pressure drag of the grooved plate is expressed as here is 45 degrees (AR = 2).  For the model description of , the experimental results of Feng et al. [45] show that the static pressure distribution in the groove can be described by an exponent or polynomial. On the basis of their work, combined with the pressure distribution in Figures  17 and 18, we assume that can be described by Equation (22).

=
where = 0, and and are adjustable variables that aim to obtain high-precision model results. Then, Equation (21) is transformed into Equation (23).
where represents the maximum pressure on the windward side, which is caused by the stagnation of the slip velocity, so can be estimated as 2 = 0.5 (24) in which is an adjustable variable and "0.5 " represents the dynamic pressure of the fluid velocity from stagnation to 0. Substituting Equations (20) and (23)  . Equation (25) shows that is related to the slip velocity and the depth of the groove. By adjusting the appropriate and values, the predicted by Equation (25) can be consistent with the numerical results. Figure 19 shows the prediction of Equation (25) with different and . Taking the fitting degree as the objective function, and -which have the largest fitting degree-are selected as the model parameters, where = 6.93 × 10 ℎ . and = 0.00955ℎ . . For the model description of p(x), the experimental results of Feng et al. [45] show that the static pressure distribution in the groove can be described by an exponent or polynomial. On the basis of their work, combined with the pressure distribution in Figures 17 and 18, we assume that p(x) can be described by Equation (22).
where C 3 = 0, and C 1 and C 2 are adjustable variables that aim to obtain high-precision model results. Then, Equation (21) is transformed into Equation (23).
where p S + 2 represents the maximum pressure on the windward side, which is caused by the stagnation of the slip velocity, so p S + 2 can be estimated as in which K p is an adjustable variable and "0.5ρ(U s ) 2 " represents the dynamic pressure of the fluid velocity from U s stagnation to 0. Substituting Equations (20) and (23) into Equation (19) yields η p as where . Equation (25) shows that η p is related to the slip velocity and the depth of the groove.
By adjusting the appropriate K 1 and K 2 values, the η p predicted by Equation (25) can be consistent with the numerical results. Figure 19 shows the prediction of Equation (25) with different K 1 and K 2 . Taking the fitting degree R 2 as the objective function, K 1 and K 2 -which have the largest fitting degree-are selected as the model parameters, where K 1 = 6.93 × 10 −6 h −1.428 and K 2 = 0.00955h −0.31646 .

Balance and to Solve the Quasi-Analytical Solution of Groove Depth ( ) in Different Inflow Velocities ( )
Simultaneous Equations (12), (18) and (25) yield Equation (26), which can predict the optimal depth and maximum depth of the transverse groove according to different inflow velocities. The steps for solving this equation are shown in Figure 20.
Step 1: We input different and to get different slip velocities of the transverse grooves.
Step 2: Input and into implicit Equation (18) to solve , and input and into Equation (25) to get .
Step 3: Taking the sum of and as the objective value. When this value is the maximum, the relationship between optimal depth and is obtained, and when this value is equal to 0, the relationship between maximum depth and is obtained. Figure 21 shows the solution when the step of U is 1 ( ∈ 1, 70 , the corresponding is calculated by Equation (2)) and the step of h is 0.0006 (ℎ ∈ 0.01, 1.2 , is calculated by corresponding h).  (12), (18) and (25) yield Equation (26), which can predict the optimal depth and maximum depth of the transverse groove according to different inflow velocities. The steps for solving this equation are shown in Figure 20.
Step 1: We input different H + and U + ∞ to get different slip velocities U + s of the transverse grooves.
Step 2: Input U + s and U + ∞ into implicit Equation (18) to solve η ν , and input U + s and H + into Equation (25) to get η p .
Step 3: Taking the sum of η ν and η p as the objective value. When this value is the maximum, the relationship between optimal depth H + opt and U + ∞ is obtained, and when this value is equal to 0, the relationship between maximum depth H + max and U + ∞ is obtained. Figure 21 shows the solution when the step of U is 1 (U ∈ [1,70], the corresponding U + ∞ is calculated by Equation (2)) and the step of h is 0.0006 (h ∈ [0.01, 1.2], H + is calculated by corresponding h).

Verification of the Quasi-Analytical Solution
In the previous section, we have established a theoretical model and solved the quas analytical solution of groove depths. In this section, firstly, the quasi-analytical solutio of this model is verified by comparing the previous research and CFD results in the pre sent with the model prediction results. Secondly, we redesign three grooved plates wit different depths based on the quasi-analytical solution and compare the drag reductio effect of the grooved plates to further verify the correctness of the model. Figure 21 shows the dimensionless depths predicted by the model versus the dimen sionless inflow velocities. The results show that the predicted optimal depths are in goo agreement with those from the present numerical simulations. However, for the predic tion of the maximum depths, the model accuracy decreases with the increase of . Tabl 5 shows the details of the prediction accuracy of this model. For the prediction of the op timal depth, the average error of the model is less than 4.35% when < 20.07 (Re 2.18 × 10 ). For the prediction of the maximum depth, the average error of the model less than 4.09% when < 19.07 (Re < 1.31 × 10 ). The main reason for the error ma be that the changes of h and Re in the CFD cases are discrete, so the presentation of th

Verification of the Quasi-Analytical Solution
In the previous section, we have established a theoretical model and solved the quasi analytical solution of groove depths. In this section, firstly, the quasi-analytical solution of this model is verified by comparing the previous research and CFD results in the pre sent with the model prediction results. Secondly, we redesign three grooved plates with different depths based on the quasi-analytical solution and compare the drag reduction effect of the grooved plates to further verify the correctness of the model. Figure 21 shows the dimensionless depths predicted by the model versus the dimen sionless inflow velocities. The results show that the predicted optimal depths are in good agreement with those from the present numerical simulations. However, for the predic tion of the maximum depths, the model accuracy decreases with the increase of . Table  5 shows the details of the prediction accuracy of this model. For the prediction of the op timal depth, the average error of the model is less than 4.35% when < 20.07 (Re < 2.18 × 10 ). For the prediction of the maximum depth, the average error of the model i less than 4.09% when < 19.07 (Re < 1.31 × 10 ). The main reason for the error may be that the changes of h and Re in the CFD cases are discrete, so the presentation of the

Verification of the Quasi-Analytical Solution
In the previous section, we have established a theoretical model and solved the quasianalytical solution of groove depths. In this section, firstly, the quasi-analytical solution of this model is verified by comparing the previous research and CFD results in the present with the model prediction results. Secondly, we redesign three grooved plates with different depths based on the quasi-analytical solution and compare the drag reduction effect of the grooved plates to further verify the correctness of the model. Figure 21 shows the dimensionless depths predicted by the model versus the dimensionless inflow velocities. The results show that the predicted optimal depths are in good agreement with those from the present numerical simulations. However, for the prediction of the maximum depths, the model accuracy decreases with the increase of U + ∞ . Table 5 shows the details of the prediction accuracy of this model. For the prediction of the optimal depth, the average error of the model is less than 4.35% when U + ∞ < 20.07 (Re < 2.18 × 10 5 ). For the prediction of the maximum depth, the average error of the model is less than 4.09% when U + ∞ < 19.07 (Re < 1.31 × 10 5 ). The main reason for the error may be that the changes of h and Re in the CFD cases are discrete, so the presentation of the optimal or maximum h in the CFD cases is not an accurate value (more details shown in Table 4).  Table 6 shows the details of the comparison between the predicted results of the model with the numerical results of Wu et al. [5] and the experimental results of Liu et al. [32]. Wu et al. [5] optimized the groove depth on the airfoil with a Reynolds number of 299,444 (U + ∞ = 20.72). The results show that when H + is 7.225, the drag reduction rate of the groove is the largest, while when H + is 21.675, the drag reduction rate of the groove approaches 0. On the premise of consistent dimensionless velocity, the prediction results of our model for the optimal and maximum depth are 8.475 and 19.290, respectively, which is close to the work of Wu et al. [5]. The inconsistency of the application object may have been the reason for the minor errors. Meanwhile, the optimal dimensionless groove depth for pipelines, H + opt = 8.488, was observed by the water tunnel experiment of Liu et al. [32] with a Reynolds number of 50,000 (U + ∞ = 17.326). In this case, the optimal dimensionless groove depth predicted by the model is 9.215, and the relative error between the model results and the experimental results is 7.8%, which qualitatively proves the correctness of the model. The inconsistency of the medium may have been the reason for the minor errors.  Figure 22, which is transformed from Figure 21, shows the relationship between the optimal and maximum groove depths and the local Reynolds numbers. Based on this model, three grooved plates have been designed at U = 4.6 m/s, as shown in Figure 23. A 160 mm smooth plate is placed in front of the grooved plate in each case. The grooved plate of Case 1 is divided into four sections, each with a length of 80 mm. According to the average local Reynolds numbers of each section, the corresponding optimal groove depths predicted by the model are 0.44, 0.35, 0.28, and 0.23 mm, respectively. On the contrary, the groove depths of Case 2 (h ≡ 0.5 mm) and Case 3 (h ≡ 0.21 mm) are calculated by the model according to the local Reynolds numbers before and after the grooved plate, respectively.      Table 7 shows that the drag reduction rate of the grooved plate (Case 1) according to the model design is the largest (12.25%) compared with other cases. On the contrary, the drag reduction rate of the grooved plate (Case 2) with the most deviation from the model design is the smallest-only 0.06%. The drag reduction rate of Case 3 is 7.83%, which is significantly less than Case 1. The reason why the drag reduction rate of Case 2 is almost zero is that the depth of the tail of the grooved plate has exceeded the predicted maximum depth (as shown in Figure 22), so these grooves at the tail cannot reduce the total drag, resulting in the decrease of the total drag reduction rate. Although the depth of the grooves in Case 3 do not meet the design of the optimal depths, it is less than the maximum groove depths predicted by the model, so Case 3 can slightly reduce the total drag.  Figure 24 compares the time-averaged entropy generation over different plates. The weaker the entropy generation, the smaller the total drag. Compared with the baseline plate, the entropy generation of the three grooved plates decreases gradually at the beginning of the groove area. From X = 0.16 to X = 0.48 in Case 1, the entropy generation at each place decreases significantly, which means that the grooves at each place play a role in reducing drag. However, in Case 2, the closer the groove is to the tail section of the grooved plate, the smaller the reduction of entropy generation, which means that the groove in the tail section may not reduce drag because the groove depth here exceeds the maximum depth shown in Figure 22. Compared with Case 2, the entropy generation of the front section (X = 0.32) of Case 3 decreases less, and the entropy generation of the tail section (X = 0.40~0.48) decreases more because the groove depth of the front section of Case 2 is closer to the optimal design depth and the groove depth of the tail section of Case 3 is closer to the optimal design depth, which further verifies the correctness of the model.

Conclusions
In this paper, we used the LES with the dynamic subgrid model to investigate the drag reduction characteristics of transverse V-grooves with different depths (h = 0.05~0.9 mm) at different Reynolds numbers ( 1.09 × 10 4 ∼ 5.44 × 10 5 ). Based on the numerical results, the physical model describing the relationship between the dimensionless depth (H + = Hu τ /υ) of the transverse groove, the dimensionless inflow velocity (U + ∞ = U ∞ /u τ ), and the drag reduction rate (η) was established to quasi-analytically solve the optimal and maximum transverse groove depth according to different Reynolds numbers. The main conclusions are summarized as follows: (1) The dimensionless groove depth (H + ) and dimensionless inflow velocity (U + ∞ ) affect the magnitude of the slip velocity (U + s ), thus driving the variation in the total viscous drag and pressure drag and thereby affecting the total drag. Therefore, the essence of balancing the benefit (the viscous drag reduction rate, η ν ) and cost (the pressure dragincrease rate, η p ) to obtain the optimal drag reduction is to obtain the optimal slip velocity by matching the depth of the groove and the inflow velocity.
(2) The relationship between U + s and η ν can be constructed by comparing the velocity profile of the grooved plate (slip) and baseline plate (no-slip), and the relationship between U + s and η p can be established by analyzing the effect of the slip velocity on the pressure distribution of the groove wall. When the value of η ν + η p reaches the maximum, the relationship between optimal depth H + opt and U + ∞ is obtained, and when this value is equal to 0, the relationship between maximum depth H + max and U + ∞ is obtained.
(3) The model results are consistent with the present numerical results and with previous data. For the solution of the optimal depth, the average error of the model is less than 4.35% when U + ∞ < 20.07 (Re < 2.18 × 10 5 ). For the solution of the maximum depth, the average error of the model is less than 4.09% when U + ∞ < 19.07 (Re <1.31 × 10 5 ).

Conclusions
In this paper, we used the LES with the dynamic subgrid model to investigate the drag reduction characteristics of transverse V-grooves with different depths (h = 0.05~0.9 mm) at different Reynolds numbers (1.09 × 10 ~5.44 × 10 ). Based on the numerical results, the physical model describing the relationship between the dimensionless depth (H = Hu υ ⁄ ) of the transverse groove, the dimensionless inflow velocity (U = U /u ), and the drag reduction rate (η) was established to quasi-analytically solve the optimal and maximum transverse groove depth according to different Reynolds numbers. The main conclusions are summarized as follows: (1) The dimensionless groove depth (H ) and dimensionless inflow velocity (U ) affect the magnitude of the slip velocity (U ), thus driving the variation in the total viscous drag and pressure drag and thereby affecting the total drag. Therefore, the essence of balancing the benefit (the viscous drag reduction rate, ) and cost (the pressure drag-increase rate, ) to obtain the optimal drag reduction is to obtain the optimal slip velocity by matching the depth of the groove and the inflow velocity.
(2) The relationship between U and can be constructed by comparing the velocity profile of the grooved plate (slip) and baseline plate (no-slip), and the relationship between and can be established by analyzing the effect of the slip velocity on the pressure distribution of the groove wall. When the value of reaches the maximum, the relationship between optimal depth and is obtained, and when this value is equal to 0, the relationship between maximum depth and is obtained. Institutional Review Board Statement: Not applicable.

Informed Consent Statement: Not applicable.
Data Availability Statement: Not applicable.