Analytic Solution of Optimal Aspect Ratio of Bionic Transverse V-Groove for Drag Reduction Based on Vorticity Kinetics

Abstract

Previous studies have implied that the AR (aspect ratio) of the transverse groove significantly affects the stability of the boundary vortex within the groove and thus drives the variation in the drag-reduction rate. However, there is no theoretical model describing the relationship between the AR and the stability of the boundary vortex, resulting in difficulty in developing a forward method to obtain the optimum AR. In this paper, the velocity potential of the groove sidewalls to the boundary vortex is innovatively described by an image vortex model, thus establishing the relationship between the AR and the induced velocity. Secondly, the velocity profile of the migration flow is obtained by decomposing the total velocity inside the groove, by which the relationship between the AR and the migration velocity is established. Finally, the analytical solution of the optimal AR ($A{R}_{opt}=2.15$) is obtained based on the kinematic condition for boundary vortex stability, i.e., the induced velocity equals the migration velocity, and the forms of boundary vortex motion at other ARs are discussed. Furthermore, the stability of the boundary vortex at the optimal AR and the corresponding optimal drag-reduction rate are verified by the large eddy simulations method. At other ARs, the motion forms of the boundary vortex are characterized by “vortex shedding” and “vortex sloshing,” respectively, and the corresponding drag-reduction rates are smaller than those for vortex stability.

1. Introduction

Reducing the friction drag of aircraft has been an active research field for several decades. An assessment showed that reducing the skin-friction drag by 1% of a Boeing 777F freighter means annual savings of around 3700 tons of kerosene [1]. Among many drag-reduction technologies, the transverse-grooved surfaces imitating dolphin skin [2,3] that are regular perpendicular to the streamwise direction have been extensively investigated for engineering applications because of their remarkable drag-reduction properties (the highest drag-reduction rate in the flat plate was up to 16% [4]) and good applicability [5,6].

There is no definite conclusion regarding the actual drag-reduction mechanism of transverse-grooved surfaces, and some different theories have been adopted. In one view, the transverse grooves on a flat wall should be classified as the ”D” and ”K” type roughness. The vortices formed within the grooves weaken the turbulence structure in the boundary layer near the wall, thus delaying the transition from laminar to turbulent flow [7,8]. The other more popular perspective indicates that the vortices (they are called boundary vortices in this paper) within transverse grooves change the sliding friction into the rolling friction at the solid–liquid interface, which is also named the “micro air-bearing phenomenon” [9,10,11,12,13,14,15]. In addition, the studies of Mariotti et al. [16,17], Pasqualetto et al. [18], Howard et al. [19], and Lang et al. [3,20] showed that the vortices formed in the transverse grooves (and small cavities) increase the momentum in the boundary layer near the wall, thus effectively controlling the flow separation. However, none of the qualitative descriptions of the drag-reduction mechanism can be applied in the design of grooves, so a large number of optimization studies on the geometric parameters of transverse grooves have to be conducted through numerical simulations and experiments in engineering applications.

Previous parametric studies have found that the shape, spacing, depth, and AR (aspect ratio, the ratio of width to height) of the transverse groove are important parameters affecting the drag-reduction rate. Cui and Fu [21] conducted a numerical study 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. Wu et al. [22,23] carried out a series of experiments and numerical simulations to explore the drag reduction induced by crescent grooves and V-grooves. The results showed that the drag-reduction rate of V-grooves is better than that of crescent grooves. The experimental and numerical results of Liu et al. [24] indicated that the drag-reduction performance in a pipeline was significantly affected by the spacing and widths of grooves, and a maximum drag-reduction rate of 3.21% could be achieved at Re = 89,500. Li et al. [25] investigated the drag-reduction characteristics of nanoporous transverse grooves through a water tunnel experiment and found that a 40% drag-reduction rate was obtained with a groove depth of 100 nm. Bhatia et al. [26] explored the effect of groove depth on drag reduction and found that the drag-reduction rate reached 10.8% when the ratio of groove depth to boundary-layer thickness was about 0.003. Liu et al. [27] proposed a new idea for numerical optimization from the perspective of entropy generation and adopted a Multi-Island Genetic Algorithm (MIGA) to optimize the geometry of grooves. These works have a certain significance for guiding engineering applications. However, conducting a large number of numerical and experimental optimization works without a theoretical basis can be time-consuming, challenging, and expensive. In order to guide the engineering design of transverse grooves, a practical alternative is to construct a theoretical model to solve the mathematical relationship between geometric parameters based on the physical conditions that affect the drag-reduction properties.

The research of Gautam [28] inspired us to solve the mathematical relationship between geometric parameters based on the stability of the boundary vortex. He conducted experiments on cavity-embedded flat plates with different ARs. The results showed that with the cavity geometry varied, the stability of the boundary vortex and the streamlines at the slip boundary also changed, which in turn affected the surface drag. In other words, it indicates that the AR affects the stability of the boundary vortex, thus driving the variation in the drag-reduction rate. Unfortunately, to the best of the authors’ knowledge, the physical relationship between the AR and the stability of the vortex within the groove has not been established. It is significant as a guide for the engineering design of the optimal AR. Therefore, the purpose of this study is to solve the optimal AR of the transverse V-groove based on the vorticity kinematics conditions of boundary vortex stability.

This paper is organized as follows. First, the influence of boundary vortex stability on drag-reduction performance is qualitatively analyzed in Section 2. Second, in Section 3, the relationship between the AR and kinematic parameters (velocity induced by image vortices and migration velocity decomposed by mainstream velocity) is established, respectively, and then the optimal AR is solved based on the equilibrium between induction velocity and migration velocity. Then, in Section 4, the vorticity kinematics characteristics of the boundary vortices in grooves with different ARs are calculated by numerical simulation, and the drag-reduction performances of grooves with different ARs are compared, which is then followed by conclusions in Section 5.

2. Influence of Boundary Vortex Stability on Drag-Reduction Performance

According to our previous research [4], the drag-reduction performance of grooves is determined by the reduction in viscous resistance and the increase in pressure drag. Moreover, we analyzed the influence of groove depth on the viscous drag-reduction rate (benefits) and the pressure drag-increase rate (costs). In this section, another key factor, the stability of the boundary vortex, that affects the compromise between benefits and costs is qualitatively analyzed to prove that it is also a necessary condition for obtaining the optimal drag-reduction rate.

Figure 1a,b show the schematic diagram of the slip surface above the transverse grooves when the boundary vortices are stable and unstable, respectively. When the boundary vortices are stable, they act as “air bearings,” separating the boundary layer from the solid wall, thus forming a smooth slip surface that results in fluid sliding over the grooved plate, while the vortices sloshing in the grooves may lead to large fluctuations in the slip surface when the boundary vortices are unstable. On the one hand, these fluctuations may cause dramatic shear between the boundary layer and the boundary vortices, and the aerodynamic configuration of the “air bearing” is destroyed, both of which are unfavorable for reducing viscous drag. On the other hand, the unsmooth slip surface aggravates the stagnation of the high-momentum fluid in the boundary layer on the groove windward side, which increases the additional pressure drag compared with the case when the boundary vortices are stable.

To further analyze the influence of boundary vortex stability on the flow within the grooves, Figure 2 compares the schematic diagrams of velocity profiles when the boundary vortex is stable and unstable. Figure 2a shows the case that the boundary vortex stabilizes on a symmetrical V-shaped transverse-grooved plate without the adverse pressure gradient and disturbance sources. In this state, the boundary vortex is stable on the centerline of the groove and acts as “fluid bearings” [9,10] on the boundary to lubricate the free stream “slip over” [9] the grooved regions. It induces a steady slip velocity (for more definitions and descriptions of slip velocity, refer to [4]), $U_{s1}$, on the horizontal surface, which reduces the velocity gradient near the wall and significantly reduces viscous resistance in comparison to a flat plate (the solid line shown in Figure 2a is the velocity profile on the centerline of the groove, and the dashed line is the velocity profile on the flat plate at the corresponding position). In addition, large-scale stagnation of fluid on the windward side can be avoided due to the “smoothness” of the slip surface, and then the increase in pressure drag is prevented [4,13,14]. When the boundary vortex is unstable, it moves upstream or downstream. Figure 2b shows the case that the boundary vortex deviates from the center of the groove and moves upstream. In this case, the slip surface fluctuates with the movement of the boundary vortex, and then the fluid stagnates over a wide range on the windward side, forming a wide high-pressure area. Moreover, assuming that the circulation remains constant whether the boundary vortex is stable or unstable, then the magnitude of slip velocity, $U_{s2}$, corresponds to the projection of $U_{s1}$ that is fictional in Figure 2b on the horizontal axis (i.e., $U_{s2}<U_{s1}$). It is detrimental to the reduction in viscous resistance because the velocity gradient near the wall is not significantly reduced. Similarly, as shown in Figure 2c, when the boundary vortex moves downstream, part of the fluid from the groove’s leeward surface moves outwards induced by the boundary vortex, forming a wide relative low-pressure area. Meanwhile, the performance of reducing the viscous resistance is also worse than that when the boundary vortex is stable due to the decrease in slip velocity ($U_{s3}<U_{s1}$). In summary, as shown in Figure 2b,c viscous resistance decreases slightly and the pressure drag increases substantially when the boundary vortex is unstable, both of which result in unsatisfactory drag-reduction performance.

The analysis above shows that maintaining the stability of the boundary vortex within grooves is the essential vorticity kinematics condition for perfect drag-reduction performance. Without considering the effects of an adverse pressure gradient, a curvature of the plate, and external perturbations, the stability of the boundary vortices is mainly related to the AR of the transverse groove [28,29]. Therefore, based on the vorticity kinematics condition of boundary vortex stability, the theoretical solution of the optimal AR is solved in Section 3.

3. Theoretical Solution of AR for Maintaining the Stability of Boundary Vortices

In the groove of a 2D plate without adverse pressure gradient and external perturbations, the stability of the boundary vortices is affected by the groove wall (which prevents the boundary vortices from moving downstream) and the migration flow (which causes the boundary vortices to move downstream) [30,31,32], as shown in Figure 3. The influence of a groove on a boundary vortex can be quantitatively described by the velocities induced by two image vortices with equal vorticity and direction opposite to that of the boundary vortex [31,32]. In Figure 3, $U_{BC}$ represents the velocity induced by image vortex B at the center of boundary vortex C, $U_{AC}$ represents the velocity induced by image vortex A at the center of boundary vortex C, and $U_m$ represents the vector sum of these two induced velocities, which is called the total induced velocity. The effect of migration flow on the boundary vortex is described by the migration velocity $U_c$.

The kinematic condition of vortex stability claims that the boundary vortex is stable when $\left|U_m\right|=U_c$ [31]. Based on this equation, the optimum AR can be derived in the procedure shown in Figure 4:

Step 1. Construct the relationship between AR and induced velocity $U_m$.

Step 2. Construct the relationship between AR and migration velocity $U_c$.

Step 3. The theoretical solution of the optimal AR is obtained based on the dynamic conditions of vortex stability ($\left|U_m\right|=U_c$).

3.1. Induced Velocity Induced by Image Vortices

To quantitatively evaluate the total induced velocity $U_m$, it is necessary to deduce the velocity field induced by the two image vortices. The cross-section of an image vortex is assumed to be a circle with a radius R surrounded by an unbounded irrotational flow, as shown in Figure 5. This vortex has a constant vorticity $\omega$ within the circle and with no vorticity outside. At point $P$, the velocity magnitude is $q$, which can be expressed as

$$q=\{\begin{array}{c}\frac{\omega r}{2}, r<R\\ \frac{\omega R^2}{2r}, r\ge R\end{array}$$

(1)

where $r$ denotes the distance between point $P$ and the center of the image vortex. Assume that the angle between the direction of flow velocity at point $P$ and the real axis of the complex plane is $\theta +\pi /2$, and then the velocity vector at point $P$ can be described as in Equation (2)

$$U+Vi=q{e}^{\left(\theta +\pi /2\right)i}$$

(2)

in which $U$ represents the velocity component in the horizontal direction, and $V$ stands for the vertical velocity component.

When the flow field outside the image vortex is assumed to be potential flow, the velocity field around the vortex can be described with a simple and general method of the complex potential function represented by $W=\phi +\psi i$, where $\phi =\phi \left(x^{\prime },y^{\prime }\right)$ denotes the potential function outside the image vortices, and $\psi =\psi \left(x^{\prime },y^{\prime }\right)$ represents the flow functions outside the image vortices. The coordinates of a point $P$ outside the image vortex are written in the form $Z=x^{\prime }+y^{\prime }i=r{e}^{\theta i}$, and then the derivative of $W$ with respect to $Z$ at point $P$ can be represented by Equation (3).

$${-\frac{dW}{dZ}|}_{Z_p}=U-Vi=q{e}^{-\left(\theta +\pi /2\right)i}$$

(3)

Because point $P$ can represent any point outside the vortex, substituting Equation (1) into Equation (3) and integrating the equation yields $W$.

$$W=\frac{{\Gamma }_{iv}i}{2\pi }lnZ$$

(4)

in which ${\Gamma }_{iv}$ represents the circulation of the image vortex, which is equal to $\omega \pi R^2$. For the vortex with circulation ${\Gamma }_{iv}$ and centered at $Z_0$, the complex potential function outside the vortex is

$$W=\frac{{\Gamma }_{iv}i}{2\pi }\ln \left(Z-Z_0\right)$$

(5)

where $Z_B=-a/2+0i$ and $Z_C=a/2+0i$ represent the coordinates of the center points of image vortices A and B, respectively, in the coordinate system shown in Figure 3. Therefore, the complex flow function of the two image vortices in a groove can be obtained by substituting the coordinates of the centers of the image vortices: $W_B=\frac{{\Gamma }_{iv}i}{2\pi }\ln \left(Z+\frac{a}{2}\right)$ $W_C=\frac{{\Gamma }_{iv}i}{2\pi }\ln \left(Z-\frac{a}{2}\right)$. In addition, the complex potential function of the flow field induced by the two image vortices is (i.e., $W_B+W_C$)

$$W_m=\frac{{\Gamma }_{iv}i}{2\pi }\ln \left(Z^2-\frac{a^2}{4}\right)$$

(6)

Now, differentiation of Equation (6) with respect to $Z$ at point $P$ ($Z_P=yi$) yields Equation (7), which represents the velocity induced by the image vortices at the centerline of the transverse grooves

$$U_{iv}=-\frac{{\Gamma }_{iv}y}{\pi \left(y^2+\frac{a^2}{4}\right)}$$

(7)

Moreover, as shown in Figure 3, the relationship between a, b, and AR can be derived as $a=\frac{4bh}{s}=\frac{4b}{AR}$, where h and s represent the groove height and width, respectively. Therefore, the total induced velocity $U_m$ at the center of the boundary vortex, i.e., the total induced velocity, is stated by Equation (8), which is determined by the circulation of the image vortex and the geometric parameters of the groove

$$U_m=-\frac{{\Gamma }_{iv}}{\pi \left(1+\frac{4}{A{R}^2}\right)b}$$

(8)

where $b\in \left(0, \frac{2s^{2}h}{s^2+4h^2}\right)$ and $AR\in \left(0,+\infty \right)$. Figure 6a,b show six representative realizations of Equations (7) and (8), respectively. When $b$ is a constant value that represents the position of the boundary vortex, the velocity profiles induced by image vortices within the groove vary with ARs, as shown in Figure 6a, which means that the velocity profile inside the groove is significantly affected by the side walls of the groove. Figure 6b shows the variation in the virtual distance $R_{vd}=\left|U_m\right|\times \pi /{\Gamma }_{iv}$ [33] that is inversely proportional to the total induced velocity with b and AR at the center of the boundary vortex. It can be observed that when AR is constant, the virtual distance $R_{vd}$ decreases with the increase in b, which means that the fuller the boundary vortex, the greater the total induced velocity.

In sum, Equation (8) establishes the relationship between the AR of the groove and the induced velocity of the image vortex, which describes the velocity field potential of the sidewall inside the groove.

3.2. Migration Velocity Decomposed by Total Velocity

From the perspective of vorticity kinematics, the total velocity profile ($U_{ms}$) inside a transverse groove is affected by the boundary vortex, image vortices, and migration flow, as described in Equation (9)

$$U_{ms}=U_{iv}+U_{bv}+U_{mv}$$

(9)

where $U_{ms}$ represents the total velocity, $U_{iv}$ is the total induced velocity that is solved by Equation (8), $U_{bv}$ denotes the velocity components induced by the boundary vortex that can be expressed as Equation (1), and $U_{mv}$ stand for the migration velocity. Therefore, the velocity profile of the migration flow within the groove can be obtained by subtracting the components of the velocity profile associated with the image vortices and boundary vortex.

Emily Jones [34] measured the total velocity profile of a transverse groove plate with AR = 1.33, h = 2 mm, and $R{e}_h=750$ (based on the groove height), as shown by the black dotted line in Figure 7a. The AR and the position of the vortex can be substituted into Equation (8) to obtain the total induced velocity of the image vortices, as shown by the red dotted line. Meanwhile, the velocity induced by the boundary vortex can be obtained according to Equation (1). Finally, the migration velocity decomposed by the total velocity in Equation (9) can be derived, as shown by the green dotted line. Similarly, the migration flow in another groove (AR = 2, h = 0.2 mm, and $R{e}_h=68$) is obtained by decomposing the corresponding total velocity profile calculated in our previous study [4], as shown in Figure 7b. It can be observed that the total velocity profiles in the groove vary with the AR, h, and $R{e}_h$, but the migration velocity profiles after decomposition are similar.

The experimental and numerical studies of Feng et al. [35] imply that both the velocity profile and pressure distribution inside the groove can be described by an exponent or polynomial function. On the basis of their work, combined with the migration velocity profiles in Figure 7, we assume that the migration velocity profile within grooves is

$$U_{mv}^{+}=C_1{\left(Y^{+}\right)}^2+C_2{Y}^{+}+C_3$$

(10)

in which $U_{mv}^{+}$ is the dimensionless form of migration velocity, and $Y^{+}$ is the dimensionless form of ordinate Y ($U_{mv}^{+}=U_{mv}/U_{\infty }$ and $Y^{+}=Y/h$, where $U_{\infty }$ denotes the mainstream velocity). $C_1$, $C_2$, and $C_3$ are parameters related to the slope of the points D and C, as shown in Figure 8. The migration velocities of points D and C refer to the velocities at the upper boundary (i.e., the slip surface) and the center of the boundary vortex, respectively.

According to the slip theory, the velocity gradient on the slip surface can be expressed as

$${\frac{\partial u}{\partial y}|}_{slip}=\frac{U_s^{+}}{L_s^{+}}$$

(11)

in which $U_s^{+}$ is the dimensionless form of slip velocity, and $L_s^{+}$ is the dimensionless form of the slip velocity length. For points D and C, the slip velocities can be regarded as $U_d^{+}$ and $U_c^{+}$, respectively, and the slip lengths can be considered as $L_d^{+}$ and $L_c^{+}$, and then the slopes $K_d$ and $K_c$ for points D and C can be expressed as

$$K_d=\frac{U_d^{+}}{L_d^{+}}, K_c=\frac{U_c^{+}}{L_c^{+}}$$

(12)

where $U_d^{+}$ and $L_d^{+}$ stand for the migration velocity of the slip surface and the distance from the slip surface to the center of the boundary vortex, respectively. Similarly, $U_c^{+}$ and $L_c^{+}$ are the migration velocity of the center of the boundary vortex and the distance from the center of the boundary vortex to the bottom of the groove, respectively. According to the geometric relationship shown in Figure 3, $L_d$ and $L_c$ can be represented by Equation (13).

$$L_d=R=\frac{1}{2}\sqrt{b^2+\frac{a^2}{4}}, L_c=\frac{b}{2}+\frac{a^2}{8b}$$

(13)

Under the established coordinate axis shown in Figure 8, Equation (14) can be obtained by substituting Equations (12) and (13), and the coordinates of point C, (0, $U_c^{+}$), into Equation (10).

$$U^{+}=\left(\frac{U_d^{+}}{2L_d^{+}{}^2}-\frac{U_c^{+}}{2L_d^{+}{L}_c^{+}}\right)\times {\left(Y^{+}\right)}^2+\frac{U_c^{+}}{L_c^{+}}{Y}^{+}+U_c^{+}$$

(14)

On the basis of Equation (14), the relationship between $U_d$ and $U_c$ can be obtained by substituting the coordinates of point D, ($L_d^{+}$,$U_d^{+}$).

$$U_c=\frac{U_d}{\left(\frac{L_d}{L_c}+2\right)}$$

(15)

In order to verify the relationship between $U_c$ and $U_d$ in Equation (15), Figure 9 shows the comparison of the experimental results and numerical calculations with the theoretical formula. It is observed that the theoretical formula is in good agreement with the results of previous studies, which further indicates that the assumption of migration velocity is reasonable.

The value of $U_d$ in Equation (15) cannot be expressed quantitatively without numerical calculation or experiment. Our previous work shows that it is related to the circulation of the boundary vortex (${\Gamma }_{bv}$) and the geometric parameters of the groove [4]. Therefore, if the boundary vortex circulation, ${\Gamma }_{bv}$, is assumed to be a known quantity, the physical relation between the ${\Gamma }_{bv}$ and $U_d$ can be established first, and then the value of $U_c$ can be measured quantitatively by ${\Gamma }_{bv}$ and geometric parameters. Moreover, the vorticity of the boundary vortex within the transverse groove can be approximated as $\omega =-\partial u/\partial y\approx -U_d/L_d$. Therefore, $U_d$ can be expressed as $U_d=-{\Gamma }_{bv}/\pi L_d$, and then substituted for Equation (15), the quantitative formula of $U_c$ can be expressed as

$$U_c=\frac{{\Gamma }_{bv}}{-\pi L_d\left(\frac{L_d}{L_c}+2\right)}$$

(16)

Equation (16) states that the physical relation between the migration velocity caused by the migration flow (i.e., $U_c$) and geometric parameters ($L_d$ and $L_c$ in Equation (16) are related to $a$ and $b$, as shown in Equation (13)).

3.3. Solution of Aspect Ratio Based on Equalization between Induced Velocity and Migration Velocity

In Section 3.2 and Section 3.3, the formulas are derived by calculating the induced velocity and migration velocity. The kinematic condition for boundary vortex stability requires that these two velocities are equal [31]; based on this equalization, the optimal AR is solved in this section.

The circulation of the image vortex is opposite to that of the boundary vortex

$${\Gamma }_{iv}=-{\Gamma }_{bv}$$

(17)

The kinematic condition that maintains the stability of the boundary vortex is

$$\left|U_m\right|=U_c$$

(18)

Substituting Equations (8), (16), and (17) into Equation (18) yields the mathematical relation between $a$ and $b$. Then, it can be obtained that the optimal AR ($A{R}_{opt}$) of the groove that maintains the stability of the boundary vortex should satisfy.

$$A{R}_{opt}=\frac{s}{h}=\frac{4b}{a}=2.15$$

(19)

It is worth noting that considering the influence of external disturbances in practical applications, the optimal AR may be a value close to this theoretical solution.

In order to further understand the influence of AR on the motion trajectory of the boundary vortex, it is necessary to discuss the cases when $A{R}_{max}>AR>A{R}_{opt}$ and $A{R}_{min}<AR<A{R}_{opt}$ from the perspective of vorticity kinematics. $A{R}_{max}$ and $A{R}_{min}$ are two critical values corresponding to the maximum AR and minimum AR, respectively. When AR exceeds these two values, the boundary vortex may not be formed inside the groove [9,29]. When $A{R}_{min}<AR<A{R}_{opt}=2.15$, Equation (18) is transformed into

$$\left|U_m\right|<U_c$$

(20)

Then, the boundary vortex moves downstream initially. It is assumed that the core of the boundary vortex moves over $\Delta x$ from point $C$ to point $C^{\prime }$, as shown in Figure 10. Therefore, the centers of the image vortices move from point $A$ and $B$ to point $A^{\prime }$ and $B^{\prime }$, respectively. The coordinates of these points are as follows:

$A\left(\frac{a}{2},0i\right)$	$A^{\prime }\left(\frac{a}{2}+\Delta x\cos 2\theta ,\Delta xi\sin 2\theta \right)$
$B\left(-\frac{a}{2},0i\right)$	$B^{\prime }\left(-\frac{a}{2}+\Delta x\cos 2\theta ,-\Delta xi\sin 2\theta \right)$
$C\left(0,bi\right)$	$C^{\prime }\left(\mathsf{\Delta }x,bi\right)$

In the representations of coordinates, $\theta$ is the angle between the windward side of the groove and the horizontal plane. Substituting the coordinates of point $A^{\prime }$, $B^{\prime }$, and $C^{\prime }$ into Equations (3)–(8) yields the complex potential function of the image vortex, which can be described as

$$W_{A\prime B\prime }=\frac{{\Gamma }_{iv}i}{2\pi }\ln [{\left(Z-\Delta x\cos 2\theta \right)}^2-{\left(\frac{a}{2}+i\Delta x\sin 2\theta \right)}^2]$$

(21)

Here, $a$ and $b$ can be replaced by Equation (22).

$$\{\begin{array}{c}a=2L_c\sin 2\theta \\ b=2L_c{\cos }^2\theta \end{array}$$

(22)

Then, the differentiation of Equation (21) with respect to Z in “$Z_{c\prime }=\Delta x+bi$” yields Equation (23)

$$-U_m+V_{m}i={\frac{dW}{dZ}|}_{Z_{c\prime }}=\frac{{\Gamma }_{iv}\left[-\Delta x\left(\cos 2\theta -1\right)i+2L_c{\cos }^2\theta \right]}{\pi \left[4\Delta x^2{\sin }^2\theta -4L_c{}^2{\cos }^2\theta \right]}$$

(23)

Therefore, $U_m$ and $V_m$ can be expressed in Equation (24)

$$\{\begin{array}{c}{U}_m=\frac{{\Gamma }_{iv}{L}_c{\cos }^2\theta }{\pi \left[2\Delta x^2{\sin }^2\theta -2L_c{}^2{\cos }^2\theta \right]}\\ V_m=-\frac{{\Gamma }_{iv}\Delta x\left(\cos 2\theta -1\right)}{\pi \left[4\Delta x^2{\sin }^2\theta -4L_c{}^2{\cos }^2\theta \right]}\end{array}$$

(24)

Two equations derived above, i.e., Equations (16) and (24), illustrate that the motion of the boundary vortex can be divided into three stages:

Stage I. Since $\Delta x\sin \theta \ll L_c\cos \theta$ and $\cos 2\theta <1$, $V<0$ in Equation (24) can be inferred. This result shows that when the boundary vortex deviates from the center of the groove and moves downstream, the vertical velocity induced by the image vortices makes the boundary vortex move toward the negative direction of the Y-axis in Figure 10b. With the boundary vortex moving to the bottom of the groove, $L_c$ decreases and $L_d$ increases in Equation (16), so dramatic decreases in $U_c$ can be inferred, which ultimately results in $\left|U_m\right|>U_c$. Therefore, the trajectory of the boundary vortex is shown as the St_1 line in Figure 10b.

Stage II. When the boundary vortex moves through the centerline of the groove, $\Delta x$ in Equation (24) is less than 0, so $V>0$ can be inferred. This result means that the vertical velocity induced by the image vortices makes the boundary vortex move toward the horizontal line. Therefore, the boundary vortex trajectory is shown as the St_2 line in Figure 10b.

Stage III. When the boundary vortex is close to the horizontal line, the effect of the migration flow on the boundary vortex is more and more prominent, and eventually the relationship between the induced velocity and the migration velocity returns to the original relationship, i.e., $\left|U_m\right|<U_c$. Therefore, the boundary vortex trajectory is shown as the St_3 line in Figure 10b.

Similarly, when $A{R}_{max}>AR>A{R}_{opt}=2.15$, Equation (18) is transformed into

$$\left|U_m\right|>U_c$$

(25)

Hence the boundary vortex will migrate upstream initially under the influence of the image vortices. It is assumed that the core of the boundary vortex moves over $\Delta x$ from point $C$ to point $C\prime$. Then, the centers of the image vortices move from points $A$ and $B$ to points $A\prime$ and $B\prime$, respectively, as shown in Figure 11. The coordinates of the centers of the boundary vortex and the image vortices are as follows:

$A\left(\frac{a}{2},0i\right)$	$A^{\prime }\left(\frac{a}{2}-\Delta x\cos 2\theta ,-\Delta xi\sin 2\theta \right)$
$B\left(-\frac{a}{2},0i\right)$	$B^{\prime }\left(-\frac{a}{2}-\Delta x\cos 2\theta ,\Delta xi\sin 2\theta \right)$
$C\left(0,bi\right)$	$C^{\prime }\left(-\mathsf{\Delta }x,bi\right)$

Consistent with the theoretical derivation above, the velocity induced by the image vortices to the boundary vortex can be expressed as

$$\{\begin{array}{c}{U}_m=\frac{{\Gamma }_{iv}{L}_c{\cos }^2\theta }{\pi \left[2\Delta x^2{\sin }^2\theta -2L_c{}^2{\cos }^2\theta \right]}\\ V_m=\frac{{\Gamma }_{iv}\Delta x\left(\cos 2\theta -1\right)}{\pi \left[4\Delta x^2{\sin }^2\theta -4L_c{}^2{\cos }^2\theta \right]}\end{array}$$

(26)

According to Equation (26), it is inferred that $V<0$ when the boundary vortex deviates from the center of the groove and moves upstream. Conversely, $V>0$ when the boundary vortex moves downstream. Therefore, the boundary vortex moves along lines St_1, St_2, and St_3 in Figure 11b. When $\theta \to 0$ (which means that the boundary vortex is difficult to form) or ${\Gamma }_{bv}\to 0$ (i.e., the boundary vortex is constantly dissipating), it can be inferred that $U_m\to 0$ and $V\to 0$, which explains why the boundary vortex initially moves upstream along the windward direction but eventually disperses with the migration flow.

When the boundary vortex moves to the left and the right limit positions, the horizontal induced velocity $U_m$ of the boundary vortex caused by the image vortices is equal to the migration velocity $U_c$ ($V\ne 0$ at this time). Then, Equation (27) is constructed from Equations (16), (24), and (26).

$$\frac{-{\Gamma }_{iv}{L}_c{\cos }^2\theta }{\pi \left[2\Delta x^2{\sin }^2\theta -2L_c{}^2{\cos }^2\theta \right]}=-\frac{{\Gamma }_{bv}}{\pi \left[\left(L_c\cos \theta -\Delta x\sin \theta \right)\left(\cos \theta +2\right)\right]}$$

(27)

When the boundary vortex fills the groove, $\left(1+\cos \theta \right)L_c\Delta AR=S$ can be deduced. Then, $\Delta x/L_c$ and $\left|\Delta x/S\right|$ can be expressed as

$$\{\begin{array}{c}\frac{\Delta x}{L_c}=\frac{{\cos }^2\theta \left(2+\cos \theta \right)-\cos \theta \sqrt{{\cos }^4\theta +4{\cos }^3\theta -4{\cos }^2\theta -16\cos \theta +16}}{4\sin \theta }\\ \left|\frac{\Delta x}{s}\right|=\left|\frac{{\cos }^2\theta \left(2+\cos \theta \right)-\cos \theta \sqrt{{\cos }^4\theta +4{\cos }^3\theta -4{\cos }^2\theta -16\cos \theta +16}}{4\sin \theta \left(1+\cos \theta \right)\Delta AR}\right|\end{array}$$

(28)

Figure 12 shows that the maximum dimensionless horizontal moving distance $\left|\Delta x\right|/s$ varies with the AR. It is worth noting that $\left|\Delta x\right|/s=0$ when the AR = 2.15, which means that the boundary vortex is stable in this case, and further proves that the above theoretical derivation is reasonable. As the AR increases from 0 to 2.15, $\left|\Delta x\right|/s$ gradually decreases to 0. When the AR increases from 2.15 to $+\infty$, $\left|\Delta x\right|/s$ gradually increases from 0 with a decreasing growth rate. With the increase in $\left|\Delta x\right|/s$, the momentum and the scale of the boundary vortex decrease. Therefore, it is reasonable to assume that it is difficult to maintain the periodic sloshing when the scale of the boundary vortex is less than $K_{1}R+K_{2}s\sin \theta$, as shown in Figure 10 and Figure 11. K₁ and K₂ are two empirical parameters equal to 0.2 and 0.1, respectively. Then, Equation (28) is derived, describing the maximum offset distance to maintain boundary vortex sloshing.

$$\left|\Delta x\right|/s\le \frac{\left(1-K_1\right)\cos \theta }{\sin \theta \left(1+\cos \theta \right)AR}-K_2$$

(29)

Combining Equations (28) and (29) yields $A{R}_{min}=0.75$ and $A{R}_{max}=6.15$, as shown in Figure 12. Considering the complex flow patterns in engineering applications, it is worth noting that the critical ARs may be the values close to theoretical solutions, as shown in the parameter band in Figure 12, which is more significant in providing a theoretical basis for groove design when considering processing errors.

In conclusion, according to the above-derived vorticity kinematics theory, without considering the adverse pressure gradient and external disturbance, the motion forms of the boundary vortex inside the groove can be divided into three forms with the variation in the AR. For AR = 2.15, i.e., region I in Figure 12, the boundary vortex is stable inside the groove. When $A{R}_{max}=6.15>AR>A{R}_{opt}=2.15$ and $A{R}_{min}=0.75<AR<A{R}_{opt}=2.15$, that is, region II in Figure 12, because the maximum horizontal moving distance is less than the limit value (blue dotted line), the boundary vortex is sloshing in the groove. When $AR<A{R}_{min}=0.75$ and $AR>A{R}_{max}=6.15$, i.e., region III in Figure 3, as the maximum horizontal movement distance is greater than the limit value, the boundary vortex is difficult to maintain with periodic sloshing inside the groove and may migrate downstream with the mainstream.

4. Numerical Verification

The numerical method with large eddy simulation (LES) was tested in the plate flow over transverse grooves in our previous studies [4]. In this section, more details and verification of this numerical method are introduced.

4.1. Numerical Methodology

The large eddy simulation (LES) with a dynamic subgrid-scale (SGS) is used in the commercial software FLUENT 18.0 to investigate the induced drag reduction and flow characteristics [36]. Moreover, a central-differencing scheme is used for spatial discretization, and a second order implicit time-stepping approach is used for temporal discretization. The space and time resolution of the numerical method is of second-order accuracy. The dimensionless physical timestep $\mathsf{\Delta }tU/h\approx 0.02$ [4] and dimensionless statistical averaging time $TU/h\approx 400$ (greater than 1000 time steps [36]) are used, where $U$ denotes the uniform velocity at the inlet and $h$ represents the depth of grooves.

The overall dimensionless size of $L_x^{+}\times L_y^{+}\times L_z^{+}=3112\times 466\times 311$ is shown schematically in Figure 13 and

Table 1. Size of the computational domain and the corresponding number of grid nodes.

		Dimensionless Parameters	Nodes
$L_x^{+}$		2490 + 311 + 311	300 + 1000 + 60
$L_y^{+}$		466	80
$L_z^{+}$		311	80
$\Delta x^{+}$	Groove	0.3	1000
	Other	<10	300 + 60
$\Delta y^{+}$		0.02~10	80
$\Delta z^{+}$		3.9	80

, where x, y, and z axes denote the streamwise, wall-normal, and spanwise directions, respectively. The size of the computational domain is larger than the minimum flow unit suggested by Jiménez and Moin [37]. The smooth walls with a dimensionless length of 2490 (160 mm) and 311 (20 mm) are located upstream and downstream of the grooved wall to prevent the propagations of pressure perturbations at the inlet and outlet, respectively. The simulated grooved wall is about 20 mm long, consisting of different symmetric V-groove profiles whose depths are 0.2 mm (the groove depths of 0.2 mm are chosen in order to ensure that the grooves have a high drag-reduction rate at both Reynolds numbers of $5.44\times {10}^4$ and $9.8\times {10}^4$) and ARs are 0.5, 1, 2, 5, and 8, where the ARs = 0.5 and 8 in region III shown in Figure 12 are selected to verify the motion form of the boundary vortex moving downstream with the mainstream, the ARs = 1 and 5 belonging to region II are used to verify the motion form of vortex sloshing, and the example of AR = 2 is used to verify the stability of the boundary vortex in region I.

The Reynolds numbers are $5.44\times {10}^4$ are $9.80\times {10}^4$, which are based on the length of the flat wall placed upstream of a grooved wall (160 mm). The no-slip condition is specified for all the solid walls, and the symmetry boundary conditions are applied to the upper and lateral sides of the computational domain [38]. At the inlet of the computational domain, an ideal gas flow with uniform velocity is set. The spectral synthesizer provides an alternative method of generating fluctuating velocity components at the inlet. In this method, fluctuating velocity components are computed by synthesizing a divergence-free velocity-vector field from the summation of Fourier harmonics (In detail, the turbulence intensity of 0.1% is set) [39,40]. Moreover, the outlet is given as the pressure outlet.

Figure 14 shows the structured mesh around the transverse V-grooves that are generated by ICEM. 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 ${\mathrm{Y}}^{+}$ is 0.02. The maximum normal grid resolution $\Delta y_{max}^{+}$ is less than 10. The streamwise grid resolution $\Delta x^{+}$ is 0.3 within the grooves, and $\Delta x_{max}^{+}$ is less than 10 in other streamwise positions. The spanwise grid resolution $\Delta z^{+}$ is 3.9 [4,36].

In order to verify that the grid resolution meets the requirements of the large eddy simulation, two streamwise grid resolutions within the groove $\Delta x^{+}$ (0.3 and 1.2) and two spanwise grid resolutions $\Delta z^{+}$ (3.9 and 5.2) are chosen to investigate the effect on the outcomes.

Table 2. Verification of grid resolution (h = 0.2 mm, $Re=5.44\times {10}^4$ ).

$\mathit{\Delta }{\mathit{x}}^{+}\left(\mathbf{Groove}\right)$	$\mathit{\Delta }{\mathit{z}}^{+}$	Drag (N)	Streamline
0.3	3.9	0.00314
	5.2	0.00313
1.2	3.9	0.00314
	5.2	0.00315

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 $\Delta x^{+}=0.3$ (groove) and $\Delta z^{+}=3.9$, which are selected for the derivation of all the other results.

The results of the grid-independent validation at five different grid-refinement levels based on the comparison of the drag-reduction rates are presented in Figure 15, with the vertical coordinates indicating the relative error in the drag-reduction rate relative to the total number of grids at 8.7 million. The results show that the relative error is only 0.087% when the number of grid cells exceeds 4.1 million. Moreover,

Table 3. Calculation of the grid convergence index by Richardson extrapolation method.

Variables	Value
$N_1$ (coarse)	2,404,420
$N_2$ (medium)	4,123,710
$N_3$ (fine)	8,703,310
$r_{12}$	1.19
$r_{23}$	1.28
p	$\approx$5.77
$g_p^{12}$	0.0860%
$g_p^{23}$	0.0159%
$GC{I}^{12}$	0.0795%
$GC{I}^{23}$	0.026%

shows the grid convergence studies with total resistance as a variable; the grid convergence indexes of $GC{I}^{12}$ and $GC{I}^{23}$ are less than 0.1%. The results mean that both the medium and fine grids meet the requirements of calculation accuracy. Given the complexity of the flow within the groove, the case of 8.7 million grid cells is used as the grid resolution to derive all other results.

Figure 16 compares the numerical calculation results with the experimental results in References [41,42]. The simulated velocity profiles over the grooved plates with a groove depth of 1.62 mm and a groove width of 3.57 mm are in good agreement with the experimental results obtained by Ahmadi-Baloutaki et al. [41] as shown in Figure 16a,b (turbulence intensity is 0.5% and 4.4%, respectively), with a maximum relative less than 3%. Moreover, as shown in Figure 16c,d, the velocity vectors obtained numerically and experimentally [42] inside the groove are identical, indicating that the CFD method can accurately capture the flow details over the grooved plate.

4.2. Stability of Boundary Vortices and Drag-Reduction Rate of Transverse Grooves with Different ARs

Figure 17 shows the coherent structures of the boundary vortices identified by ${\lambda }_{ci}$ in three dimensions. It can be observed that despite the random perturbations given in the computational domain (see Figure 17a), the boundary vortices can still be identified in the grooves (see Figure 17b). Since the transverse groove is isotropic in the spanwise direction, and the coherent structures of the boundary vortices have no obvious change in the spanwise direction and the unsteady phenomenon within the groove mainly changes in the streamwise direction, all the 3D simulation results are shown for the cross-sections located at 50% of the spanwise of the domain, as shown in Figure 17b.

Figure 18, Figure 19 and Figure 20 show the evolution of the instantaneous spanwise vorticity in the grooves with different ARs at Re = 5.44E4, where the ARs = 0.5 and 8 (Figure 18) at region III, the ARs = 1 and 5 (Figure 19) belonging to region II, and the AR = 2 (Figure 20) at region I were predicted theoretically as shown in Figure 12.

When the AR = 0.5, there are two vortices with different scales in the groove, as shown in Figure 18a. The small-scale vortex is shed from the shear layer and moves along the windward side to the bottom of the groove. The scale of the vortex becomes larger because it is stretched during the clockwise shaking in the groove and eventually dissipates in the groove. Similarly, such a vortex shedding phenomenon also occurs in the groove with AR = 8, as shown in Figure 18b. After shedding from the shear layer at the leeward side, the vortex gradually migrates downstream with the mainstream and then gradually dissipates near the windward side. Since the flow phenomena in the grooves with AR = 0.5 and AR = 8 are qualitatively consistent with the shedding of the vortex street, it is classified as the “vortex shedding” phenomenon. Because this vortex shedding is unsteady, it is difficult to completely describe using the theory of quasi-steady vorticity kinematics. However, this phenomenon justifies the range of region III predicted by our theory. In this range, the vortex periodically moves downstream due to an initial $\Delta x$ after the vortex shedding is greater than the critical value and finally dissipates near the windward side.

Figure 19a,b show the periodic sloshing of the boundary vortex in the groove with AR = 1 and AR = 5, respectively. This periodic development process can be divided into two stages. In the first stage, the boundary vortex moves along the windward side to the groove bottom, corresponding to t = 0.039 s to t = 0.043 s in Figure 19a and t = 0.035 to t = 0.037 in Figure 19b. In the second stage, the boundary vortex moves clockwise along the leeward side to near the horizontal line, corresponding to t = 0.044 s to t = 0.050 s in Figure 19a and t = 0.040 to t = 0.045 in Figure 19b. The periodic motion of the boundary vortex in the groove with AR = 1 and AR = 5 is classified as the “vortex sloshing” phenomenon.

In contrast, when the AR is 2, the boundary vortex stays stable in the groove, as shown in Figure 20. In this case, the slip surface separating the boundary layer from the solid wall is smooth, so the fluid flows smoothly through the groove, resulting in less viscous loss caused by the friction between the fluid and the solid compared with the baseline plate and less momentum loss caused by vortex shedding or sloshing compared with other grooved plates.

In order to further distinguish the three motion forms of the boundary vortex, the time histories of the dimensionless vertical velocities ($V_s^{+}$) at the intersection point of the groove centerline and the horizontal line at $Re=5.44\times {10}^4$ and $Re=9.8\times {10}^4$ are compared in Figure 21. It is observed that the dimensionless vertical velocities at this point fluctuate periodically, which implies that the vortex motion is periodic and continuous. Moreover, the variance in vertical velocity is compared in Figure 22 to explore the variation in the $V_s^{+}$ fluctuation degree with Reynolds number. The results show that the variances of AR = 0.5 and AR = 8 are close, which means that the influence of boundary vortex shedding in these ARs on velocity fluctuations near the slip surface is similar. Meanwhile, the variances in dimensionless vertical velocities are close for AR = 1 and AR = 5, which further suggests that the motion form of the boundary vortices at these ARs is similar. Furthermore, the fluctuation degree of $V_s^{+}$ is always minimal at AR = 2, implying that the boundary vortices stay stable in the grooves at $Re=5.44\times {10}^4$ and $Re=9.8\times {10}^4$.

Figure 23a illustrates the variation in the total drag-reduction rate with the AR at different Reynolds numbers. The total drag-reduction rate is defined as

$$\eta =\frac{F_G-F_R}{F_R}$$

(30)

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 in AR at $Re=5.44\times {10}^4$ and $Re=9.8\times {10}^4$. The maximum drag-reduction rate of 11.59% appears at AR = 2, and the minimum drag-reduction rate appears at AR = 0.5 and 8. In particular, the drag-reduction rate approaches 0 at AR = 8. This further suggests that the AR directly affects the stability of the boundary vortex, which in turn drives the variation in the groove drag-reduction rate. As mentioned in Section 3.3, the optimum value of the AR may vary slightly given the complex flow phenomena in engineering applications, but it is a value close to 2.15.

The total drag of the grooved plate consists of viscous drag ($F_{GV}$) and pressure drag ($F_{GP}$), which are expressed by Equations (31) and (32), respectively. $F_{GV}$ and $F_{GP}$ are determined by calculating the corresponding local stress, namely, shear ($\mathsf{\tau }$) and pressure (P) at the wall, and integrating the projected stress in the drag direction ($\overline{n}$, that is X direction) along the wetted wall ($l_s$).

$$F_{GV}={{\displaystyle \int }}_0^{l_s}{\left[\overline{n}·\tau \right]}_s·xdx$$

(31)

$$F_{GP}={{\displaystyle \int }}_0^{l_s}{\left[\overline{n}·P\right]}_s·xdx$$

(32)

Therefore, the total drag-reduction rate is transformed into Equation (33)

$$\eta =\frac{F_{GV}-F_R}{F_R}+\frac{F_{GP}}{F_R}={\eta }_{\nu }+{\eta }_p$$

(33)

Here, “${\eta }_{\nu }=\left(F_{GV}-F_R\right)/F_R$” denotes the reduction rate of viscous drag, and “${\eta }_p=F_{GP}/F_R$” denotes the increased rate of pressure drag. Figure 23b,c show the variation in ${\eta }_{\nu }$ and ${\eta }_p$ with the AR, respectively. The results show that the absolute value of ${\eta }_{\nu }$ increases first and then decreases, and ${\eta }_p$ is minimal at AR = 2, which means that the increase in additional pressure drag is minimal when the boundary vortex is stable. These results are consistent with the analysis in Section 2.

5. Conclusions

The relationship between the AR of the transverse groove and the motion of boundary vortices is established based on vortex kinematics, and the analytical solution of the optimal AR for maintaining the stability of the boundary vortices is solved. Moreover, the optimal AR is validated by large eddy simulations (LES), and the motion of the boundary vortex for other ARs is analyzed. The main conclusions are as follows.

(1)

The velocity potential of the groove sidewalls to the boundary vortex is described by an image vortex model, thus establishing the relationship between the AR and the induced velocity. Secondly, the velocity profile of the migration flow is obtained by decomposing the total velocity inside the groove, by which the relationship between the AR and the migration velocity is established. Finally, the analytical solution of the optimal AR ($A{R}_{opt}=2.15$) is obtained based on the kinetic conditions (i.e., the induced velocity is equal to the migration velocity) of the boundary vortex stability and the value of the critical ARs ($A{R}_{min}=0.75$ and $A{R}_{max}=6.15$) for which the boundary vortex can slosh inside the groove is obtained. Without considering the adverse pressure gradient and external disturbance, the motion forms of the boundary vortex inside the groove can be divided into three forms with the variation in the AR.

(2)

The theoretical model for solving the optimal AR ($A{R}_{opt}$) and critical ARs ($A{R}_{min}$ and $A{R}_{max}$) is validated by investigating the motion of the boundary vortices and the drag-reduction rate of the groove for ARs of 0.5, 1, 2, 5, and 8 with large eddy simulations. For AR = 2, the boundary vortex is stable inside the groove, corresponding to the maximum drag-reduction rate. When the AR is closer to 2, i.e., AR = 1 and AR = 5 (corresponds to the interval $A{R}_{min}<AR<A{R}_{opt}$ and $AR<A{R}_{opt}<A{R}_{max}$), the boundary vortices slosh periodically inside the groove and the magnitude of the vertical velocity fluctuations is similar in both cases. This periodic motion of the boundary vortex in the groove is classified as the “vortex sloshing” phenomenon. When the AR is far from 2, i.e., AR = 0.5 and AR = 8 (corresponds to the interval $AR<A{R}_{min}$ and $AR>A{R}_{max}$), the boundary vortices are shed from the shear layer at the leeward side and migrate downstream with the mainstream, which is classified as the “vortex shedding” phenomenon and corresponds to the minimum drag-reduction rate.

It is worth emphasizing that, given the complex flow phenomena in engineering applications, there may be a slight uncertainty in the actual critical and optimal AR, but both should be close to the theoretical solution. Therefore, the following experimental work will focus further on this possible difference.