Shear Stress in Turbulent Flow over Smooth Surfaces

Avraham Rami Oron

Independent Researcher, Givatayim, Israel.

DOI: 10.4236/oalib.1113015                

Abstract

This paper presents a general equation for shear stress in turbulent flow over smooth surfaces at Reynolds numbers ranging from 2 × 10³ to 6.5 × 10⁵ (based on the boundary layer thickness and the velocity at its edge). The equation is based on a mathematical analysis using the experimental turbulent flow data in a smooth circular pipe as a reference. The general shear stress equation presented here was tested by solving the momentum equation for several surfaces, including a circular pipe, rectangular ducts, flat plates with and without suction, a very long cylinder perpendicular to the stream, and supersonic turbulent flow over a flat plate and along a circular cylinder. It showed good compatibility with the experimental data.

Keywords

Turbulent, Shear

1. Introduction

Turbulence occurs when the stream’s kinetic energy is too high, and the viscosity cannot dissipate all of it into heat.

In the presence of a wall, all stress effects occur within a boundary layer, where the shear stress is the most dominant and the rest are negligible. The shear stress itself is divided into laminar and turbulent flows. Additionally, as Reynold’s number increases, the effect of the laminar shear stress decreases until it becomes negligible.

The research in turbulent flow computation began when J. Boussinesq [1] formulated the eddy viscosity theorem. In analogy with the coefficient of laminar flow, he suggested a turbulent viscosity that depends on the velocity derivative components of disturbance in the x and y directions. Prandtl [2] proposed the mixing length theory, in which the shearing stress depends on the square of the first derivative of the main velocity, multiplied by the square of the mixing length and the density. This model yields good results, but each case’s mixing length is specific and needs to be obtained experimentally.

Based on the empirical data of flow in a pipe, Schlichting [3] suggested calculating the flow field by the integral of the momentum equation. He assumed that inside the boundary layer, the velocity profile and the shearing stress are the same as those of a pipe, where the thickness of the boundary layer replaces the radius of the pipe. He inserted these equations into the momentum equation, and its integration yielded the dependence of the boundary layer thickness on the distance from the leading edge x, and the shear stress on the surface.

The most updated model is the k-w (Wilcox [4]). It predicts turbulence using two differential equations—the first (k) for the turbulence kinetic energy and the second (w) for the specific dissipation rate of the turbulent kinetic energy into internal thermal energy.

In this paper, we derive the shear stress equation for validation, examine turbulent flow over a smooth surface, and find exact solutions for the momentum equation. Subsequently, we compare these solutions to experimental data of turbulent flow.

2. Derivation of the Shear Stress Equation

2.1. Modeling the Framework

The fundamental of obtaining the shear stress of turbulent flow over a smooth surface is to assume a universal function relating it to the first derivative of the main velocity and apply it to known experimental shear stress data. A mathematical analysis yields the unique possibility for a general equation. This equation is used to solve the momentum equation in many cases. The equation can be assumed to be general if the various solutions are compatible with the experimental data.

The experimental case chosen is the flow in a circular pipe, which was profoundly investigated. The shearing stress has been used to calculate the shear stress in many other smooth surfaces.

2.2. Flow through a Smooth Circular Pipe—Experimental Data

H. Schlichting [3] summarized the turbulent flow in a circular pipe, a topic that many researchers have explored. The empirical results for this case are:

$\frac{u}{U}={\left(\frac{y}{r}\right)}^{\frac{1}{n}}$ (2.1)

$\frac{U_m}{U}=\frac{n^2}{\left(n+1\right)\cdot \left(n+0.5\right)}$ (2.2)

${\tau }_w=\frac{ϵ}{8}\cdot \rho \cdot U_m^2$ (2.3)

$\frac{1}{\sqrt{ϵ}}=2\cdot \text{log}\left(\frac{U_m r}{\nu }\cdot \sqrt{ϵ}\right)-0.2$ (2.4)

where r is the radius of the pipe, $U_m$ is the average velocity (flow rate per pipe area), U is the velocity in the axis of the pipe, n is the exponent depending on $\frac{U_{m}r}{\nu }$ , ${\tau }_w$ is the shear stress on the surface, and ϵ is the resistance coefficient.

Table 1 presents the relation of n to $\frac{U_{m}r}{\nu }$ and the conversion of the dependence of the shear stress to $\frac{Ur}{\nu }$ .

Eqs. (2.1-2.4) yield the experimental shearing stress. A good approximation of it is:

$\frac{{\tau }_w}{\rho \cdot U^2}=0.01425\cdot {\left(\frac{\nu }{U\cdot r}\right)}^{\frac{1}{5}}$ (2.5)

A comparison of the experimental shearing stress to that obtained from Eq. (2.5) is provided in Table 1.

Table 1. Shear stress in eq. (1.5) vs. experimental data for various $\frac{U_{m}r}{\nu }$ values.

$\frac{U_{m}r}{\nu }$	2∙103	1.15∙104	5.5∙104	5.5∙105	1∙106
n	6	6.6	7	8.8	10
$\frac{Ur}{\nu }$	2.53∙103	1.43∙104	6.73∙104	6.47∙105	1.16∙106
$\frac{{\tau }_w}{\rho U^2}$ data	0.00312	0.00204	0.00147	0.00103	0.000975
$\frac{{\tau }_w}{\rho U^2}$ eq.(1.5)	0.00297	0.00210	0.00154	0.00098	0.000874

Also, the shearing stress along the vertical direction (perpendicular to the pipe axis) is:

$\tau ={\tau }_w\left(1-\frac{y}{r}\right)$ (2.6)

2.3. Mathematical Analysis

Under the assumption that there is a general equation for the shearing stress of turbulent flow over smooth surfaces, we can define:

$\mu \frac{\partial u}{\partial y}=\tau \cdot \frac{{\partial }^{2}Q}{\partial y^2}$ (2.7)

where x and y are the coordinates in the directions along and perpendicular to the surface, µ is the viscosity, and Q is a valid general function that depends only on x and y.

Since Q is a general function, eq. (2.7) is also valid in the case of turbulent flow in a pipe.

Thus, in this case, using eq. (2.6) yields:

$\mu \frac{\partial u}{\partial y}={\tau }_w\cdot \left(1-\frac{y}{r}\right)\cdot \frac{{\partial }^{2}Q}{\partial y^2}$ (2.8)

Integrating eq. (2.8) from y = 0 to y = r yields:

This equation gives:

$Q\left(r,x\right)=\frac{1}{0.01425}{\left(\frac{\mu }{\rho \cdot U}\right)}^{4/5}\cdot r^{6/5}$ (2.10)

The only possibility that Q would be a general function is if:

$Q\left(y,x\right)=\frac{1}{0.01425}{\left(\frac{\mu }{\rho \cdot U}\right)}^{4/5}\cdot y^{6/5}$ (2.11)

From this, we obtain:

$\frac{{\partial }^{2}Q}{\partial y^2}=\frac{1}{0.0594}{\left(\frac{\mu }{\rho Uy}\right)}^{4/5}$ (2.12)

Eq. (2.12), together with eq. (2.7) yields an equation for the shearing stress of incompressible flow (see eq. (3.1)). However, to make it valid also in the case of compressible flow, it will be:

$\tau =0.0594{\left(\frac{U{\displaystyle {\int }_0^y\rho \text{d}y}}{\mu }\right)}^{4/5}\cdot \mu \frac{\partial u}{\partial y}$ (2.13)

It should be noted that since eq. (2.13) is an integral of the momentum equation, it cannot detect the thin laminar sublayer near the wall.

2.4. Principle of Separated Flow Fields

The principle of separated flow fields enables the solution of 3D cases. It states that in the case of a flow field over a segmented surface, it is divided into sub-flow fields, so each segment has its own flow field, which depends on the distance from the segment and the undisturbed stream, i.e., $u=u\left(y,U\right)$ . Every two close fields are separated by a thin separation zone so that the velocity on both sides is equal and the velocity gradients become equal inside it. This principle has been applied to the flow of non-circular pipes and is supported by empirical data.

3. Exact Solutions of the Momentum Equation

3.1. Incompressible Turbulent Flow in Smooth Circular Pipes

In the case of turbulence in a pipe, the shear stress in the y direction is given by comparison of this equation with eq. (3.1) yields:

$\tau =\lambda \cdot \rho \cdot {\left(\frac{Uy}{\nu }\right)}^{4/5}\nu \frac{\partial u}{\partial y}$ $\left(\lambda =0.0594\right)$ (3.1)

Comparison of this equation with eq. (2.6) yields:

$\lambda \rho {\left(\frac{Uy}{\nu }\right)}^{4/5}\nu \frac{\partial u}{\partial y}={\tau }_w\left(1-\frac{y}{r}\right)$ (3.2)

The solution of eq. (3.2) is:

$\frac{u}{U}=\frac{1}{5}{\left(\frac{y}{r}\right)}^{1/5}\left(6-\frac{y}{r}\right)$ (3.3)

A comparison of the velocity in eq. (3.3) to that in eq. (2.1) is presented in Table 2 and exhibits small differences.

Table 2. Calculated vs. experimental velocities for various $\frac{y}{r}$ values.

$\frac{y}{r}$	0	0.2	0.4	0.6	0.8	1.0
$\frac{u}{U}$ eq. (3.3)	0	0.841	0.932	0.975	0.995	1.000
$\frac{u}{U}-{\left(\frac{y}{r}\right)}^{\frac{1}{7}}$	0	0.046	0.055	0.045	0.026	0.000
$\frac{u}{U}-{\left(\frac{y}{r}\right)}^{\frac{1}{9}}$	0	0.005	0.029	0.030	0.019	0.000

For the shear stress on the surface of the pipe, we obtain:

$\frac{{\tau }_w}{\rho U^2}=0.01425{\left(\frac{\nu }{U\cdot r}\right)}^{1/5}$ (3.4)

This equation is identical to eq. (2.5).

3.2. Incompressible Turbulent Flow in Rectangular Ducts

Turbulent flow in noncircular ducts has been investigated in numerous studies (Schlichting, [3]). The experiments yield, among others, two important phenomena:

• The curves of constant velocity are, principally, parallel to the nearest wall. This led to the above-mentioned assumption of the separated flow fields.

• The shear stress on the surface around the cross section is uniform. Thus, the constant velocity curve is also a constant shear stress.

A scheme of a rectangular duct with wall lengths of 2 H and 2 h so that $H\ge h$ is described in Figure 1.

Figure 1. Rectangular duct cross-section.

According to the principle of separated flow fields, two pairs of flow fields—#1 and #2-are separated by 4 + 1 separation zones. The 4 zones extend on the cross angles between the walls. Thus, we obtain the following:

The distance from the walls to the constant velocity curve is the same, i.e.:

$y_1=y_2=y$ (3.5)

The shear stress and the pressure gradient are related by:

$\tau \cdot \left(H+h-2y\right)=-\frac{\text{Δ}P}{\text{Δ}x}\left(h-y\right)\left(H-y\right)$ (3.6)

Eq. (3.6) can be divided into 2 equations:

${\tau }_w=-\frac{\Delta P}{\Delta x}\frac{Hh}{\left(H+h\right)}$ (3.7)

And

$\tau ={\tau }_w\frac{\left(H-y\right)\left(h-y\right)}{Hh}\frac{\left(H+h\right)}{\left(H+h-2y\right)}=\lambda {\left(\frac{Uy}{\nu }\right)}^{\frac{4}{5}}\rho \nu \frac{\text{d}u}{\text{d}y}$ (3.8)

A numerical solution of eq.(3.8) yields, for each ratio of H/h:

$\frac{u}{U}=\frac{1}{5}{\left(\frac{y}{h}\right)}^{1/5}\left(6-\frac{y}{h}\right)$ (3.9)

And

$\frac{{\tau }_w}{\rho U^2}=0.01425{\left(\frac{\nu }{U\cdot h}\right)}^{1/5}$ (3.10)

These equations are compatible with the experimental law of the hydraulic pipe.

4. Incompressible Turbulent Flow in the Boundary Layer

4.1. Governing Equations

The momentum equation in this case is:

$u\frac{\partial u}{\partial x}+v\frac{\partial u}{\partial y}=U\frac{dU}{dx}+\lambda \frac{\partial }{\partial y}\left[{\left(\frac{Uy}{\gamma }\right)}^{4/5}\nu \frac{\partial u}{\partial y}\right]$ (4.1)

while the continuity equation yields:

$u=\frac{\partial \phi }{\partial y}$ $v=-\frac{\partial \phi }{\partial x}$ (4.2)

Assuming that the velocity profile is semi-similar in shape, we can write:

$\phi =U\delta F\left(\eta ,X\right)$ (4.3)

where:

$\eta ={\left(\frac{y}{\delta }\right)}^{1/5}$ and $X=x$ (4.4)

δ: Typical length in the y direction

Setting:

$f=\frac{u}{U}$ (4.5)

So that:

$F={\displaystyle {\int }_0^{\eta }5\cdot {\eta }^4\cdot f\cdot \text{d}\eta }$ (4.6)

Inserting eqs. (4.2 – 4.6) into eq. (4.1) and separating the variables yields 2 equations. The first one, which defines the length δ, is:

$\lambda \frac{U}{\nu }={\left(\frac{U\cdot \delta }{\nu }\right)}^{1/5}\frac{\text{d}}{\text{d}x}\left(\frac{U\cdot \delta }{\nu }\right)$ (4.7)

After integration:

$\frac{U\delta }{\nu }=0.111{\left(\frac{1}{\nu }\cdot {\displaystyle {\int }_0^{x}U\text{d}X}\right)}^{5/6}$ (4.8)

The second one, which describes the relative velocity f, is:

$\begin{array}{l}\frac{{\partial }^{2}f}{\partial {\eta }^2}+5F\frac{\partial f}{\partial \eta }+25{\eta }^4\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}\left(1-f^2\right)\\ =5\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \left(5{\eta }^{4}f\frac{\partial f}{\partial X}-\frac{\partial f}{\partial \eta }\frac{\partial F}{\partial X}\right)\end{array}$ (4.9)

Here, the prime ' denotes differentiation for X.

The boundary condition for eq. (4.9) are:

$f\left(0,X\right)=0$ and $f\left(\infty ,X\right)=1$ (4.10)

After these transformations, the shear stress is:

$\tau =\frac{\lambda }{5}\rho U^2{\left(\frac{\nu }{U\delta }\right)}^{1/5}\frac{\partial f}{\partial \eta }$ (4.11)

4.2. Flow over a Flat Plate Parallel to the Stream

In this case, U is constant.

Inserting it into eq. (4.8) yields:

$\frac{\delta }{x}=0.111{\left(\frac{\nu }{Ux}\right)}^{\frac{1}{6}}$ (4.12)

And since $U^{\prime }=0$ , we obtain the turbulent flow over the flat plate equation:

$\frac{{\text{d}}^{2}f}{\text{d}{\eta }^2}+5\frac{\text{d}f}{\text{d}\eta }F=0$ (4.13)

Eq. (3.23) will be solved by re-integration, i.e., we assume an initial $f_0$ by a function that fulfills as many boundary conditions as possible at $\eta =0$ and $\eta =\infty$ . When inserting it into the equation, we obtain f₁ and so on until we acquire 2 functions, f_n, and f_n+1, that are close enough.

Since $f=\alpha \eta -\beta {\eta }^8+\cdots$ , the first function is:

$f_0={\displaystyle {\int }_0^{\eta }\left(a+b\cdot {\eta }^7\right)\exp \left(-k{\eta }^7\right)\text{d}\eta }$ (4.14)

b is calculated by the boundary condition at $\eta =\infty$ :

$b=\frac{1-{\displaystyle {\int }_0^{\eta }a\cdot \exp \left(-k{\eta }^7\right)\text{d}\eta }}{{\displaystyle {\int }_0^{\eta }{\eta }^7\exp \left(-k{\eta }^7\right)\text{d}\eta }}$ (4.15)

The next step is

$\frac{{\text{d}}^2{f}_1}{\text{d}{\eta }^2}=-5\frac{\text{d}{f}_0}{\text{d}\eta }{F}_0$ (4.16)

Integration of eq. (4.16) yields

$\frac{\text{d}{f}_1}{\text{d}\eta }=a-25\cdot \left[f_0\cdot {\displaystyle {\int }_0^{\eta }{f}_0\cdot {\eta }^4\text{d}\eta }-{\displaystyle {\int }_0^{\eta }{f}_0^2\cdot {\eta }^4\text{d}\eta }\right]$ (4.17)

And an additional integration gives $f_1$ .

The boundary condition of $f_1$ and $\frac{\text{d}{f}_1}{\text{d}\eta }$ at $\eta =\infty$ yields

$a=0.984,b=0\text{and}k=0.56$ (4.18)

The comparison of $f_1$ to $f_0$ is presented in Table 3.

Table 3. $f_0$ compared to $f_1$ vs. $\eta$ .

$\eta$	0	0.2	0.4	0.6	0.8	1.0	1.2	1.4
$f_0$	0	0.197	0.394	0.589	0.776	0.924	0.991	1
$f_1$	0	0.197	0.395	0.591	0.777	0.925	0.991	1

As can be seen, the deviation of $f_1$ from $f_0$ is less than 0.2%. Thus, the relative velocity is:

$f=0.984{\displaystyle {\int }_0^{\eta }\exp \left(-0.56{\eta }^7\right)\text{d}\eta }$ (4.19)

And

$\frac{\text{d}f\left(0\right)}{\text{d}\eta }=0.984$ (4.20)

The calculated f was compared to the empirical one, which follows the power rule $\frac{u}{U}={\left(\frac{y}{{\delta }_u}\right)}^{1/7}$ with $\frac{{\delta }_u}{x}=0.37{\left(\frac{\nu }{U\cdot x}\right)}^{1/5}$ at a range up to

$\frac{U\cdot x}{\nu }=6.4\times {10}^7$ .

This comparison is shown in Figure 2.

$\frac{{\tau }_w}{\rho U^2}=\frac{\lambda }{{\left(1.2\cdot \lambda \right)}^{1/6}}{\left(\frac{\nu }{Ux}\right)}^{1/6}\frac{0.984}{5}=0.0182{\left(\frac{\nu }{Ux}\right)}^{1/6}$ (4.21)

Table 4 shows the difference between eq. (4.21) and the experimental data (Schlichting [3]).

Figure 2. Calculated relative velocity f compared to the empirical one. The shearing stress on the surface is calculated by inserting eq. (4.12) and eq. (3.30) into eq. (3.11).

Table 4. Shear stress derived from eq. (4.21) versus experimental data for various 𝑈∞𝑥𝜈 values.

$\frac{U\cdot x}{\nu }$	5∙105	1∙106	5∙106	1∙107	5∙107
$\frac{{\tau }_w}{\rho U^2}$ data	0.00215	0.00187	0.00135	0.00118	0.00085
$\frac{{\tau }_w}{\rho U^2}$ eq. (3.32)	0.00204	0.00182	0.00139	0.00124	0.00095

4.3. Flow over Two Vertical Flat Plates (Both Parallel to the Stream)

We assume two vertical flat plates #1 and #2, as described in Figure 3.

Figure 3. The representation of flow over two vertical flat plates.

The flow fields inside the boundary layers are

$f_i=0.984{\displaystyle {\int }_0^{{\eta }_i}\exp \left(-0.56{\eta }^7\right)\text{d}\eta }$ (4.22)

where

$i=1,2$ (4.23)

The separation zone extends where the velocity of the two fields is equal, thus:

$\frac{y_1}{{\delta }_1}=\frac{y_2}{{\delta }_2}$ (4.24)

The angle between flat plate 1 and the separation $\theta$ zone is:

$\text{tang}\left(\theta \right)=\frac{y_1}{y_2}=\frac{{\delta }_1}{{\delta }_2}={\left(\frac{x_1}{x_2}\right)}^{5/6}$ (4.25)

As can be seen, $\theta =90˚$ is the leading edge of flat plate 2, and as $x_2$ increases, $\theta$ decreases until it stabilizes on the cross angle between the two flat plates.

4.4. Flow over Wedges

The potential flow over a wedge is

$U=U1\cdot X^{\frac{\theta }{\pi -\theta }}$ (4.26)

where $\theta$ is half of the wedge angle (see Figure 4)

Figure 4. Presentation of flow over a wedge.

Eqs. (4.8) and (4.26) yield

$\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}=1.2\cdot \frac{\theta }{\pi }$ (4.27)

The momentum equation is

$\frac{{\text{d}}^{2}f}{\text{d}{\eta }^2}+25\cdot \frac{\text{d}f}{\text{d}\eta }\cdot {\displaystyle {\int }_0^{\eta }{\eta }^{4}f\text{d}\eta }+25\cdot \left(1.2\frac{\theta }{\pi }\right)\cdot {\eta }^4\cdot \left(1-f^2\right)=0$ (4.28)

Eq. (4.28) will be solved using the same method used in the flat plate case. However, since $\frac{{\text{d}}^{2}f\left(0\right)}{\text{d}{\eta }^2}=-25\cdot \left(1.2\frac{\theta }{\pi }\right)\cdot {\eta }^4$ . $f$ is assumed to be

$f={\displaystyle {\int }_0^{\eta }\left[a-5\cdot \left(1.2\frac{\vartheta }{\pi }\right)\cdot {\eta }^5+b\cdot {\eta }^7\right]\exp \left(-k\cdot {\eta }^7\right)\text{d}\eta }$ (4.29)

where

$b=\frac{1-{\displaystyle {\int }_0^{\infty }\left[a-5\cdot \left(1.2\frac{\vartheta }{\pi }\right)\cdot {\eta }^5\right]\exp \left(-k\cdot {\eta }^7\right)\text{d}\eta }}{{\displaystyle {\int }_0^{\infty }{\eta }^7\cdot \exp \left(-k\cdot {\eta }^7\right)\text{d}\eta }}$ (4.30)

Then, calculating $\frac{\text{d}{f}_1}{\text{d}\eta }$ by

$\frac{\text{d}{f}_1}{\text{d}\eta }=a-25\cdot \left[f\cdot {\displaystyle {\int }_0^{\eta }{\eta }^{4}f\text{d}\eta }-{\displaystyle {\int }_0^{\eta }{\eta }^4{f}^2\text{d}\eta }+1.2\frac{\theta }{\pi }\cdot {\displaystyle {\int }_0^{\eta }{\eta }^4\cdot \left(1-f^2\right)\text{d}\eta }\right]$ (4.31)

Additional integration of eq.(4.31) gives $f_1$ .

The boundary conditions at $\eta =\infty$ that are $\frac{\text{d}{f}_1}{\text{d}\eta }=0$ and $f_1=1$ , give k, a, and b.

Table 5 presents k, a, and b for each $\theta$ .

Table 5. k, a, and b vs. $\theta$ .

$\theta$ deg.	15	30	45	60	75	90
$\frac{U\delta }{{\left(U\delta \right)}^{\prime }}\frac{U^{\prime }}{U}$	0.1	0.2	0.3	0.4	0.5	0.6
k	0.64	0.81	0.98	1.15	1.33	1.50
a	1.044	1.092	1.131	1.165	1.194	1.220
b	0.441	0.803	1.059	1.247	1.394	1.516

Figure 5 and Figure 6 present $f$ and $f_1$ for $\theta =30˚$ and $60˚$

Figure 5. $f_0$ vs. $f_1$ for $\theta =30˚$ .

Figure 6. $f_0$ vs. $f_1$ for $\theta =60˚$ .

5. Approximate Solutions in 2D Incompressible Turbulent Flow

5.1. Method of the Solution

The integral method can be used to obtain an approximate solution. This method yields a quick approximation solution without a step-by-step calculation along the vertical axis. Instead, we assume a function for the relative velocity f that fulfills as many boundary conditions as possible on the surface and at the end of the boundary layer, and integrate the momentum equation along the vertical axis.

The boundary conditions of the relative velocity f are:

$\eta =0;f=0;\frac{{\partial }^{2}f}{\partial {\eta }^2}=-25{\eta }^4\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}$ (5.1)

$\eta =\infty ;f=1;\frac{{\partial }^{n}f}{\partial {\eta }^n}=0\text{for}n\ge 1$ (5.2)

To fulfill the boundary conditions, together with the results of flow over a flat plate, the following profile for the relative velocity is assumed:

$f={\displaystyle {\int }_0^{\eta }\left(a-5{\eta }^5\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}+b{\eta }^7\right)\exp \left(-k\cdot {\eta }^7\right)\text{d}\eta }$ (5.3)

The solution of eq.(4.28) yields the value of k at the range

$0.6\ge \frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}\ge -0.27$

(with a deviation of ∓5%).

$k=\{\begin{array}{ll}0.56+0.15\cdot \frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}\hfill & \text{if}0\ge \frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}\hfill \\ 0.56+0.8\cdot \frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}\hfill & \text{if}0.1\ge \frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}\ge 0\hfill \\ 0.468+1.72\cdot \frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}\hfill & \text{otherwise}\hfill \end{array}$ (5.4)

a and b are calculated by the boundary condition at $\eta =\infty$ :

$b=\frac{1-{\displaystyle {\int }_0^{\infty }\left(a-5{\eta }^5\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}\right)\exp \left(-k\cdot {\eta }^7\right)\text{d}\eta }}{{\displaystyle {\int }_0^{\infty }{\eta }^7\exp \left(-k\cdot {\eta }^7\right)\text{d}\eta }}$ (5.5)

And by the integral of the momentum eq. (3.20):

$25\left[K+\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}S+\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{\text{d}K}{\text{d}x}\right]-a=0$ (5.6)

where

$K={\displaystyle {\int }_0^{\infty }\left(f-f^2\right){\eta }^4\text{d}\eta }$ (5.7)

And

$S={\displaystyle {\int }_0^{\infty }\left(1-f^2\right){\eta }^4\text{d}\eta }$ (5.8)

To find the initial condition, the leading edge of the surface is assumed to be a flat plate tangent to the surface, as described in Figure 7.

Figure 7. A description of the leading edge of the surface.

Thus, the leading edge is calculated as a bending flat plate, i.e.:

$\frac{\text{d}K\left(0\right)}{\text{d}X}=0$ (5.9)

The derivative of the function $K$ at point $X_{i+1}$ can be calculated as:

${\left(\frac{\text{d}K}{\text{d}X}\right)}_{i+1}=2\frac{K_{i+1}-K_i}{X_{i+1}-X_i}-{\left(\frac{\text{d}K}{\text{d}X}\right)}_i$ (5.10)

5.2. Turbulent Flow over a Very Long Cylinder

A scheme of the system is shown in Figure 8.

Figure 8. A presentation of flow over a long cylinder.

In the case of flow over a long cylinder, the velocity $U$ on the edge of the boundary layer is given by:

$U=2U_{\infty }\sin \left(\theta \right)$ (5.11)

where $U_{\infty }$ is the velocity far ahead of the cylinder and $\theta =X/r$ .

Eq. (4.8) yields:

$\frac{\delta }{r}=0.111{\left(\frac{\nu }{2U_{\infty }r}\right)}^{1/5}\cdot \frac{{\left(1-\cos \theta \right)}^{5/6}}{\sin \theta }$ (5.12)

And eq. (4.8) together with eq. (5.12)

$\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}=1.2\frac{\cos \left(\theta \right)}{1+\cos \left(\theta \right)}$ (5.13)

Thus, the relative velocity in the boundary layer is:

$f={\displaystyle {\int }_0^{\eta }\left(a-6{\eta }^5\frac{\cos \left(\theta \right)}{1+\cos \left(\theta \right)}+b{\eta }^7\right)\exp \left(-k\cdot {\eta }^7\right)\text{d}\eta }$ (5.14)

Inserting eq. (5.15) in eqs. (5.7-5.8) enables to calculate eq. (5.6) numerically (${\theta }_{i+1}-{\theta }_i=\pi /36$ ). The results are presented in Table 6.

Table 6. a vs. the angle of location θ (degrees).

θ (deg)	0	30	45	60	75	90	100
a	1.22	1.21	1.20	1.18	1.14	1.04	0.94
θ (deg)	105	110	115	120	125	128.2
a	0.87	0.79	0.67	0.52	0.30	0

Figure 9 shows the relative velocity $f\left(\eta \right)$ and its first derivative which presents the shear stress at an angle of 128.2 deg.

Figure 9. The relative velocity $f\left(\eta \right)$ and $\frac{\partial f}{\partial \eta }$ at an angle of 128.2˚.

The separation angle, 128.2˚, is similar to the separation angle that was found in an experimental study by Willy Z. Sadeh and Daniel Sharon [5] published in the NASA contractor report 3622. Their results are shown in Figure 10.

Figure 10. NASA report on the separation angle of turbulent flow over a cylinder (Sadeh, W.Z., Saharon, [5]).

5.3. Shear Stress of Turbulent Flow over a Slim Aerofoil

The parameter $\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}$ is the most influential on the shear stress. Figure 11 shows the shear stress represented by the first derivative of the relative velocity a in the cases of flow over a cylinder and flow over wedges vs. $\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}$ at ranges 0.6 to −0.2.

Figure 11. The parameter a in cases of flow over a cylinder and wedges vs. $\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\cdot \frac{U^{\prime }}{U}$ .

Assuming that the aerofoil’s leading edge is cylindrical and that $\frac{\left(U\delta \right)}{{\left(U\delta \right)}^{\prime }}\frac{U^{\prime }}{U}$ is decreased monotonically up to −5˚ at the rear edge, we get the shearing stress with a maximum deviation of 4%.

6. Supersonic Turbulent Flow of Air at Zero Pressure Gradient

This case includes all the cases where the flow consists of a constant main flow and a secondary flow, which is insignificant for the shearing stress. Thus, the shearing stress can be calculated as a flow at zero pressure gradient.

The flow equations will be solved under the following assumptions:

a) The total Prandtl number is a unit (Pt = 1).

b) The wall is adiabatic.

The viscosity is linear with the temperatures $\mu ={\mu }_w\frac{T}{T_w}$ .

The momentum equation is

$\rho u\frac{\partial u}{\partial x}+\rho v\frac{\partial u}{\partial y}=\frac{\partial \tau }{\partial y}$ (6.1)

While the continuous equation for 2-dimensional flow yields

$\rho u=\frac{\partial \phi }{\partial y}$ $\rho v=-\frac{\partial \phi }{\partial x}$ (6.2)

Under assumptions a and b, the energy equation is

$Cp{T}_w=CpT+\frac{1}{2}{u}^2$ (6.3)

Also, the equation of state yields for a perfect gas at constant pressure

$\rho T=\frac{P}{R}$ (6.4)

Setting now

${\eta }^5={\displaystyle {\int }_0^y\frac{\rho }{{\rho }_w}\frac{\text{d}y}{\delta }}$ (6.5)

And

${\eta }^5={\displaystyle {\int }_0^y\frac{\rho }{{\rho }_w}\frac{\text{d}y}{\delta }}$ (6.6)

So that

$\phi ={\displaystyle {\int }_0^y\rho u\text{d}y}={\rho }_{w}U\delta {\displaystyle {\int }_0^{\eta }f\cdot 5{\eta }^4\text{d}\eta }$ (6.7)

And the energy equation (6.3)

$\frac{T_w}{T_{\infty }}=1+0.2M^2$ (6.8)

yields

$\frac{T}{T_w}=1-\left(1-\frac{T_{\infty }}{T_w}\right)f^2=1-\frac{M^2}{5+M^2}{f}^2$ (6.9)

Inserting these expressions into the momentum equation (5.1) yields

$-\left(\frac{U{\delta }^{\prime }}{{\nu }_w}\right)\frac{\text{d}f}{\text{d}\eta }{\displaystyle {\int }_0^{\eta }f5{\eta }^4\text{d}\eta }=\frac{\partial }{\partial \eta }\left[{\left(\frac{U\delta }{{\nu }_w}\right)}^{\frac{4}{5}}\frac{\lambda }{5\delta }\frac{T_w}{T}{\left(\frac{T}{I_w}\right)}^{\frac{1}{5}}\frac{\text{d}f}{\text{d}\eta }\right]$ (6.10)

Eq. (6.10) can be divided into 2 equations—one depends on x only and the second on η only.

The first one is

$\frac{U{\delta }^{\prime }}{{\nu }_w}=\lambda {\left(\frac{U\delta }{{\nu }_w}\right)}^{\frac{4}{5}}\cdot \frac{1}{\delta }$ (6.11)

And after integration

$\frac{\delta }{x}={\left(1.2\lambda \right)}^{\frac{5}{6}}{\left(\frac{{\nu }_w}{Ux}\right)}^{\frac{1}{6}}=0.111{\left(\frac{{\nu }_w}{Ux}\right)}^{\frac{1}{6}}$ (6.12)

The second equation is

$\frac{\text{d}}{\text{d}\eta }\left[\frac{\frac{\text{d}f}{\text{d}\eta }}{{\left(1-\frac{M^2}{5+M^2}{f}^2\right)}^{0.8}}\right]=-\frac{\text{d}f}{\text{d}\eta }\cdot {\displaystyle {\int }_0^{\eta }25f{\eta }^4\text{d}\eta }$ (6.13)

The boundary conditions of eq. (6.13) are $f\left(0\right)=0$ and $f\left(\infty \right)=1$ .

Integration of eq. (6.13) yields

$\frac{\frac{\text{d}f}{\text{d}\eta }}{{\left(1-\frac{M^2}{5+M^2}{f}^2\right)}^{0.8}}=a-25f{\displaystyle {\int }_0^{\eta }{\eta }^{4}f\text{d}\eta }+25{\displaystyle {\int }_0^{\eta }{\eta }^4{f}^2\text{d}\eta }$ (6.14)

And an additional integration of the equation. (6.14) yields

${\displaystyle {\int }_0^f\frac{\text{d}f}{{\left(1-\frac{M^2}{5+M^2}{f}^2\right)}^{0.8}}}={\displaystyle {\int }_0^{\eta }\left[a-25f{\displaystyle {\int }_0^{\eta }{\eta }^{4}f\text{d}\eta }+25{\displaystyle {\int }_0^{\eta }{\eta }^2{f}^2\text{d}\eta }\right]\text{d}\eta }$ (6.15)

Eq. (6.13) can be solved numerically, but a more practical solution is to assume

$f={\displaystyle {\int }_0^{\eta }a\exp \left(-k\cdot {\eta }^7-m\cdot {\eta }^2\right)\text{d}\eta }$ (6.16)

The values of a, k, and m are calculated by equations (6.14), (6.15), and (6.16) at $\eta =\infty$ as follows: Eqs. (6.16) and (6.14) yield

$a=\frac{1}{{\displaystyle {\int }_0^{\infty }a\exp \left(-k\cdot {\eta }^7-m\cdot {\eta }^2\right)\text{d}\eta }}$ (6.17)

$a-25{\displaystyle {\int }_0^{\infty }{\eta }^2\cdot \left(f-f^2\right)\text{d}\eta }=0$ (6.18)

Then the left side of eq. (6.15) for $f=1$ is equal to the right side for $\eta =\infty$ .

The values of a, k, and b for some Mach numbers are presented in Table 7.

Table 7. The values of a, k, and m for some Mach numbers.

M	0	0.5	1	1.5	2	2.5	3	3.5
a	0.984	1.003	1.023	1.050	1.096	1.135	1.172	1.207
k	0.56	0.51	0.46	0.40	0.31	0.245	0.192	1.15
m	0	0.085	0.1754	0.291	0.48	0.634	0.775	0.9

The compatibility of f according to eq. (6.16) to that of eq. (6.15) is presented in Table 8 for M = 3.

Table 8. f per eq.(6.16) vs. f per eq.(6.15) for M = 3.

η	0.2	0.4	0.6	0.8	1.0	1.2	3
f per eq, (6.16)	0.232	0.450	0.643	0.800	0.914	0.978	1
f per eq. (6.15)	0.230	0.452	0.647	0.807	0.9205	0.9805	1.001

The distance from the surface, y, is calculated by

$\frac{y}{\delta }={\displaystyle {\int }_0^{\eta }5{\eta }^4\left(1-\frac{M^2}{5+M^2}{f}^{2\eta }\right)\text{d}\eta }$ (6.19)

The shearing stress on the surface is

${\tau }_w=\frac{0.0594\cdot a}{5}{\rho }_w{U}^2{\left(\frac{{\nu }_w}{Ux}\right)}^{0.2}$ (6.20)

6.1. 4.5 Mach Turbulent Flow over a Flat Plate

The previous equations were compared to experimental data of “4.5 Mach turbulent flow over a flat plate” published by NPARC-Alliance Validation Archive [6].

The data of this flow are as follows:

Mach number: 4.512.

Pressure: 0.0974 psia (671.6 Pa).

$T_{\infty }$ : 108.8˚R (60.5˚K).

X (distance of the measurement): 1.79 ft. (0.546 m).

The thickness of the boundary layer, ${\delta }_u$ , and the velocity profile are presented in Figure 12 and Figure 13.

Figure 12. Thickness of the boundary layer.

Figure 13. Velocity profile vs. distance from the surface.

Based on this empirical data, the following values are calculated

$T_w=T_{\infty }\left(1+M^2\right)=307˚$

${\rho }_w=\frac{P}{R{T}_w}=0.00765\text{kg}/{\text{m}}^{\text{3}}$

${\mu }_w=1.82\times {10}^{-5}\text{kg}/\text{m}\cdot \text{s}$ (from air properties calculator)

${\nu }_w=\frac{{\mu }_w}{{\rho }_w}=2.34\times {10}^{-3}{\text{m}}^{\text{2}}/\text{s}$

$U=M\cdot \sqrt{\gamma R{T}_{\infty }}=702\text{m}/\text{s}$

$\delta =0.111x{\left(\frac{{\nu }_w}{Ux}\right)}^{\frac{1}{6}}=0.0082\text{m}=8.2\text{mm}$

And Figure 7 and Figure 8 yield, by measuring ${\delta }_u=9\text{mm}$ .

Under these values, a, k, and m were found as follows:

a = 1.259.

k = 0.1.

m = 1.077.

Thus, the relative velocity is

$f=1.259{\displaystyle {\int }_0^{\eta }\exp \left(-0.1\cdot {\eta }^7-1.077\cdot {\eta }^2\right)\text{d}\eta }$ (6.21)

The calculated relative velocity f was compared to the empirical one that followsthe power profile $\frac{u}{U}={\left(\frac{y}{{\delta }_u}\right)}^{1/7}$ . This comparison is shown in Figure 14.

Figure 14. Calculated relative velocity vs. empirical one (for M4.5).

6.2. 2.4 Mach Turbulent Flow along a Circular Cylinder

Robert M. O’Donnell [7] investigated experimentally at Mach number 2.41, among the rest, the turbulent boundary layer thickness and the velocity profile. The experiments were carried out along the exterior surface of a circular cylinder. The range of the Reynolds number was $5.8\times {10}^6\ge \frac{U_{\infty }\cdot x}{{\rho }_{\infty }}\ge 0.3\times {10}^6$ . He found that the velocity profile is

$\frac{{\delta }_2}{x}=0.026\cdot {\left(\frac{{\rho }_{\infty }}{U\cdot x}\right)}^{1/5}$ (6.23)

And by converting from ${\delta }_2$ to ${\delta }_u$

$\frac{{\delta }_u}{x}=0.375\cdot {\left(\frac{{\rho }_{\infty }}{U\cdot x}\right)}^{1/5}$ (6.24)

And

$\frac{u}{U}={\left(\frac{y}{{\delta }_u}\right)}^{1/7}$ (6.25)

Calculating now the flow parameters for 2.41 Mach number yields the relative velocity

$f=1.126{\displaystyle {\int }_0^{\eta }\exp \left(-0.26\cdot {\eta }^7-0.5995\cdot {\eta }^2\right)\text{d}\eta }$ (6.26)

The typical length in the y direction $\delta$ is given by

$\frac{\delta }{x}=0.111{\left(\frac{{\nu }_w}{Ux}\right)}^{1/6}=0.144{\left(\frac{{\nu }_{\infty }}{Ux}\right)}^{1/6}$ (6.27)

The distance from the surface is given by

$\frac{y}{{\delta }_u}=\left(\frac{y}{\delta }\right)\left(\frac{\delta }{{\delta }_u}\right)=0.384\left(\frac{y}{\delta }\right){\left(\frac{Ux}{{\nu }_{\infty }}\right)}^{1/30}$ (6.28)

Setting eq. (6.28) into eq. (6.25) gives the experimental velocity.

The calculated relative velocity is compared to the experimental one in Figure 15.

Figure 15. Calculated relative velocity vs. empirical one. (for M2.4)

6.3. Supersonic Turbulent Flow on a Cylindrical Cone

A scheme of the system is shown in Figure 16.

Figure 16. Presentation of flow on a cylindrical cone.

According to the inviscid supersonic flow theory on a cylindrical cone, the velocity along the surface is constant. i.e. $\frac{\text{d}U}{\text{d}x}=0$ . Thus, the momentum and the energy that are written in eq. (6.1) and (6.3) are valid in the present case. However, since the cone is a body of revolution where the radius r is linear with the x coordinate, the continuous equation is:

$\frac{\partial }{\partial x}\left(x\rho u\right)+\frac{\partial }{\partial x}\left(x\rho v\right)=0$ (6.29)

Setting now

${\eta }^5={\displaystyle {\int }_0^{\eta }\frac{\rho }{{\rho }_w}\frac{\text{d}y}{\delta }}$ (6.30)

$\frac{u}{U}=f\left(\eta \right)$ (6.31)

And insert them, together with eq. (6.9). into the momentum equation, we get

$-\left(\frac{U{\left(x\delta \right)}^{\prime }}{x{\nu }_w}\right)\frac{\text{d}f}{\text{d}\eta }{\displaystyle {\int }_0^{\eta }f5{\eta }^4\text{d}\eta }=\frac{\partial }{\partial \eta }\left[{\left(\frac{U\delta }{{\nu }_w}\right)}^{\frac{4}{5}}\frac{\lambda }{5\delta }\frac{T_w}{T}{\left(\frac{T}{I_w}\right)}^{\frac{1}{5}}\frac{\text{d}f}{\text{d}\eta }\right]$ (6.32)

Eq. (6.32) can be divided into 2 equations, one depends on x only and the second on η only.

The first one is

$\frac{\text{d}}{\text{d}x}\left(\frac{Ux\delta }{{\nu }_w}\right)=\lambda {\left(\frac{U\delta }{{\nu }_w}\right)}^{\frac{4}{5}}\cdot \frac{x}{\delta }$ (6.33)

After integration,

$\frac{\delta }{x}=0.0574{\left(\frac{{\nu }_w}{Ux}\right)}^{\frac{1}{6}}$ (6.34)

The second equation is

$\frac{\text{d}}{\text{d}\eta }\left[\frac{\frac{\text{d}f}{\text{d}\eta }}{{\left(1-\frac{M^2}{5+M^2}{f}^2\right)}^{0.8}}\right]=-\frac{\text{d}f}{\text{d}\eta }\cdot {\displaystyle {\int }_0^{\eta }25f{\eta }^4\text{d}\eta }$ (6.35)

Which is identified as eq. (6.13) and is solved in the same.

7. Conclusions

This article presents and tests a general shear stress equation for turbulent flow over smooth surfaces. In some cases, the equation yields valid results, flow in a circular pipe, in rectangular ducts, over a flat plate parallel to the stream, the separation angle of flow over a circular cylinder, supersonic flow of air on a flat plate, and along a cylinder and turbulent flow with suction in low and high velocity.

Although these results are encouraging, many more cases must be tested until this model is approved.

Conflicts of Interest

The author declares no conflicts of interest.

List of Symbols

x	coordinate direction along the surface
y	coordinate perpendicular to the surface
$\lambda$	coefficient, equal to 0.0594
$\rho$	density
U	velocity on the edge of the boundary layer in the x direction
u	velocity inside the boundary layer in the x direction
f	relative velocity - u/U
v	velocity in the y direction
$\tau$	total shear stress
$\mu$	viscosity
$\nu$	kinematic viscosity
$\delta$	typical length in the y direction
${\delta }_u$	boundary layer thickness
${\delta }_2$	momentum thickness (${\displaystyle {\int }_0^{\infty }\frac{u}{U}\cdot \left(1-\frac{u}{U}\right)\text{d}y}$ )
P	absolute pressure
R	gas constant (285.7 kJ/kg 1˚C for air)
Cp	specific heat at constant pressure
Cv	specific heat in constant volume
$\gamma$	Cp/Cv (1.4 for air)
M	Mach number on the edge of the boundary layer
Subscript
w	refers to the surface
∞	refers to the flow far away from the surface

Conflicts of Interest

The author declares no conflicts of interest.

References

[1]	Boussinesq, J. (1877) Théorie de l’Ecoulement Tourbillant. Memoirs of the American Academy of Arts and Sciences, 26, 46.
[2]	Prandtl, L. (1926) Ueber Die Ausgebildete Turbulenz. Proceedings 2nd International Congress Applied Mechanics, Zurich, 12-17 September 1926, 62.
[3]	Schlichting, H. (1968) Boundary-Layer Theory. 6th Edition, Springer.[CrossRef]
[4]	Wilcox, D.C. (2008) Formulation of the K-W Turbulence Model Revisited. AIAA Journal, 46, 2823-2838.[CrossRef]
[5]	Sadeh, W.Z. and Saharon, D.B. (1982) Turbulence Effect on Crossflow around a Circular Cylinder at Subcritical Reynolds Numbers.
[6]	NASA (2012) Mach 4.5 Turbulent Flow over a Flat Plate. NPARC Alliance Validation Archive.
[7]	O’Donnell, R.M. (1954) Experimental Investigation at a Mach Number of 2.41 of Average Skin-Friction Coefficients and Velocity Profiles for Laminar and Turbulent Boundary Layers and an Assessment of Probe Effects. NACA Technical Note 3122.