The Solely Use of the “Hydraulic-Diameter” Concept Is Not Sufficient to Describe Correctly Non-Symmetric Flow in Conduits ()

Christos Krimizis-Tsatsoulis^{}

Alchimica S.A. R&D department, Athens, Greece.

**DOI: **10.4236/ojfd.2021.114013
PDF
HTML XML
172
Downloads
996
Views
Citations

Alchimica S.A. R&D department, Athens, Greece.

A well-known cornerstone in fluid mechanics is the equations that relate the friction factor to the Reynolds number obtained from the measurements in cylindrical cross-sectional tubes. The extension of these equations to different geometries failed to give reliable results. The introduction of the Hydraulic Diameter has fixed this issue particularly for the square ducts. However, for non-symmetric flows, as in concentric annuli, the discrepancies were unacceptable. Several attempts have been made to fix these problems with finally the introduction of a new concept like, “Laminar Equivalent Hydraulic Diameter” or “Efficient Hydraulic Diameter” provided satisfactory results. This approach seems to have fixed the problem and hence has been widely accepted. Nevertheless, it is based on a non-robust theoretical argument. In the present paper, it has been demonstrated that the solely use of the “Hydraulic Diameter” concept is insufficient to describe non-symmetric flows as in concentric annuli. It appears the need to use the Z axis component of the skew driving force for the laminar flow and the parameter λ for the turbulent one. At the same time, instead, it has been shown that in the case of flow in square and rectangular ducts, the “Hydraulic Diameter” is sufficient to describe it. In this case, the flow is practically symmetric. Moreover, several new straightforward equations are provided, which simplify a lot dealing with non-cylindrical cross-sectional conduits. In doing so, the concept of “Eigenvectors-Eigenvalues” has been implemented. This theoretical approach could help to simplify other non-symmetric cases in fluid dynamics. To mention, “Flow past immersed non-symmetric bodies”, “Flow in curved conduits” etc.

Keywords

Fluid Mechanics, Hydraulic Diameter, Flow in Annuli, Flow in Ducts, Flow of Two Adjacent Immiscible Fluids, Eigenvectors-Eigenvalues

Share and Cite:

Krimizis-Tsatsoulis, C. (2021) The Solely Use of the “Hydraulic-Diameter” Concept Is Not Sufficient to Describe Correctly Non-Symmetric Flow in Conduits. *Open Journal of Fluid Dynamics*, **11**, 210-250. doi: 10.4236/ojfd.2021.114013.

1. Introduction

There are many milestones in the fluid dynamics history. To mention the masterpiece of Navier-Stokes’ equation and the extensive work by many to understand the turbulent flow and provide the engineers with reliable design methods.

Among these are the well-known equations provided by Prandtl in 1935 for smooth pipes and Colebrook for rough pipes [1] [2]

$\frac{1}{\sqrt{f}}=2.0\mathrm{log}\left(R{e}_{{d}_{h}}\sqrt{f}\right)-0.8$

$\frac{1}{\sqrt{f}}=-2.0\mathrm{log}\left(\frac{\epsilon /d}{3.7}+\frac{2.51}{R{e}_{{d}_{h}}\sqrt{f}}\right)$ (1)

Given *Re*, these implicit equations to the friction factor *f* offer a way to calculate the pressure drop in a cylindrical tube under laminar and/or turbulent flow.

Although Equation (1) is accurate for cylindrical tubes, all attempts to extend its use to pipes and channels of different geometry have failed to provide satisfactory accuracy.

For example, the use of Equation (1) for concentric annuli shows a scattering from −25% to +35%. Attempts to cure this problem even by introducing the concept of the “hydraulic diameter” did not gave consistent results particularly for low Reynold’s’ numbers [3] [4] [5] [6].

As a matter of fact, this situation was to be expected for the following reasons:

First, Reynolds’ number is defined as a unique dimensionless combination of the governing parameters of the flow, *i*.*e*., if two or more developed flows share identical Reynolds’ numbers, then the final velocity distribution should be identical.

However, the definition of *Re* includes a linear quantity directly dependent on the geometry of the conduit. For instance, for a flow in a cylindrical pipe, the characteristic linear quantity is undoubtedly the diameter of the pipe. Hence, geometric similarity is guaranteed for all cylindrical pipes. The same reasoning applies to square ducts, where the characteristic linear dimension is just the side.

On the contrary, for other geometrical sections, like concentric annuli or rectangular ducts, a more specific definition is needed to keep the geometric similarity, specifically the *K* parameter for the annuli and the ratio *a*/*b* for the sides of a rectangular duct.

Therefore, for geometries which are other than cylindrical, the use of Equation (1) is wrong. Nevertheless, the engineers were reluctant to drop the use of Equation (1) and miss the benefit from the available work. They tried hard to attenuate the problem by introducing the idea of the “Hydraulic Diameter”. The results were acceptable for the square ducts, but unacceptable for non-symmetrical flows. Furthermore, it is unrealistic to produce an equation similar to Equation (1) for any geometric shape. Hence, the problem to relate successfully Equation (1) with other geometric shapes remains open.

On the same track, the work of Jones O. C. Jr. and Leung J. C. M [3] [4] introduced the idea of the “Laminar Equivalent Diameter”. This method is supposed to give excellent correlation to Equation (1) and so has been widely accepted [2].

In the present paper, it will be demonstrated that the main reason of these discrepancies is the non-symmetry of the flow (*i*.*e*., the skewness of the velocity distribution) rather than the lack of geometric similarity. In fact, the geometric similarity is cured well enough by just using the concept of the Hydraulic Diameter.

2. The Method

In the present paper, the approach to the flow in non-circular conduits problem is completely theoretical and based on the classical theoretical mechanics. The approach is different and the use of the Navier-Stokes’ equation and the relative concepts of Momentum Transfer, have been deliberately avoided. Finally, the concept of “Eigenvectors-Eigenvalues” has been implemented. The idea behind this decision is that by using a completely different route than the standard one the final comparison of the results will be fair.

First, a stress tensor is defined, and the relative Eigenvectors-Eigenvalues have been calculated. Consequently, the stress tensor has been diagonalized to find the direction in the fluid where shear stresses are nullified. What is left is just driving forces. Next, an equation is proposed, which provides the average velocity *V* (m/sec) as a function of the pressure gradient
$\nabla p$ (N/m^{3}). Hence, for a given average velocity the pressure gradient is calculated, or vice-versa. The results are compared with the ones provided in the literature, which make use of the widely accepted “Laminar Equivalent Diameter” concept [4]. For the laminar flow the exact analytical solutions, based on Navier stokes equation, are used for comparison. Moreover, it is assumed that the results of the work [3] [4] based on experimental data are corelating well enough with Equation (1). Hence, the results of the present work are compared with the results in [3] [4].

In the following illustrative examples this original approach becomes clearer.

2.1. Flow through an Annulus-Laminar Flow

As a first attempt, the case of flow through a concentric annulus (presented in Figure 1) is analyzed. In this case “non-symmetry” prevails.

Consider the case, $K=0.3$, ${R}_{o}=0.05\text{\hspace{0.17em}}\text{m}$, $\nabla p=0.1\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$ and water as the fluid $\mu ={10}^{-3}\text{kg}/\text{m}\cdot \mathrm{sec}$ and $\rho ={10}^{3}\text{kg}/{\text{m}}^{\text{3}}$

Calculating,
$\lambda $. *From now on we nominate **R _{o}*

Given the parameter *K*, the dimensionless parameter
$\lambda $ depends only to the geometry of the annulus and nothing more. It is calculated by the following formula,

Figure 1. Flow in a concentric annulus.

${\lambda}^{2}=\frac{1-{K}^{2}}{\mathrm{ln}\frac{1}{{K}^{2}}}$ (2)

So, $\lambda =0.615$.

Equation (2)* practically defines
${\lambda}^{2}$ as the logarithmic mean between the quantities* (*R*^{2})*and *(*K*^{2}*R*^{2}) *i*.*e*.,* the logarithmic mean of the first moment for the two wetted areas* (see also Appendix A).

Now there is need to calculate the stress on the walls ${\tau}^{w}$. Attention is required in this case as there are two different ${\tau}^{w}$. Hence, the inner wall stress as ${\tau}_{i}^{w}$ and the outer one as ${\tau}_{o}^{w}$ are designated.

Given that on the virtual cylindrical surface of radius $\lambda R$, the shear stresses are zero (maximum velocity), it is allowed to make an independent balance to the forces.

$2\pi RL{\tau}_{o}^{w}\propto \pi \left({R}^{2}-{\lambda}^{2}{R}^{2}\right)\Delta p$ (3)

$2\pi KRL{\tau}_{i}^{w}\propto \pi \left({\lambda}^{2}{R}^{2}-{K}^{2}{R}^{2}\right)\Delta p$ (4)

By dividing Equation (3) with Equation (4), the ratio of the two wall stresses is provided,

$\frac{{\tau}_{o}^{w}}{{\tau}_{i}^{w}}=K\frac{1-{\lambda}^{2}}{{\lambda}^{2}-{K}^{2}}$ (5)

At the same time, by applying a total force balance,

$2\pi RL{\tau}_{o}^{w}+2\pi KRL{\tau}_{i}^{w}=\Delta p\pi \left({R}^{2}-{K}^{2}{R}^{2}\right)$ $\iff $

${\tau}_{o}^{w}+K{\tau}_{i}^{w}=\frac{\nabla p}{2}R\left(1-{K}^{2}\right)=0.002275\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$ (6)

From Equation (5) and Equation (6), the wall stresses for this particular case are,

${\tau}_{i}^{w}=0.0024\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$ and*
${\tau}_{o}^{w}=0.00155\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$ *

As expected, ${\tau}_{i}^{w}>{\tau}_{o}^{w}$ since inner tubular surface is smaller.

The aim is to calculate the average velocity in this annulus. If the flow is a fully developed turbulent one, almost every point in the core of the fluid should share the same average velocity. In other words, in the case of the fully developed turbulent flow, the average velocity is not just a mathematical quantity, but has a physical meaning as well.

On the contrary, in the case of laminar flow, the average velocity is just a mathematical quantity, however, useful. In the core of the fluid there is a velocity distribution. In fact, due to the velocity variation there are shear stresses in the fluid. Nevertheless, a mass flow rate will be achieved at the end directly related to a theoretical average velocity.

Following the average velocity’s *V* definition as an invariant quantity in the core of the whole fluid, there is another quantity as well invariant for incompressible fluids.

$\frac{1}{2}\rho {V}^{2}=\text{KineticEnergy}$

At the same time, any element in the liquid’s core experiences stresses, which can be organized into a tensor form. There are tensile and shear stresses, which in this case, do not deform, just move the element.

$\sigma =\left|\begin{array}{cc}{\sigma}_{1}& {\tau}_{12}\\ {\tau}_{21}& {\sigma}_{2}\end{array}\right|$

Considering the element is touching the walls where there are the maximum shear stresses,

$\sigma =\left|\begin{array}{cc}-p& {\tau}_{i}^{w}\\ {\tau}_{o}^{w}& -p\end{array}\right|$

Moreover, it is understandable that the absolute magnitude of the pressure at a certain section (e.g., $Z=0$ ) Inside the tube, it has absolutely no influence on the flow behavior for an incompressible fluid. Hence, this pressure tensor can be subtracted,

$\sigma =\left|\begin{array}{cc}-p& {\tau}_{i}^{w}\\ {\tau}_{o}^{w}& -p\end{array}\right|-\left|\begin{array}{cc}-p& 0\\ 0& -p\end{array}\right|=\left|\begin{array}{cc}0& {\tau}_{i}^{w}\\ {\tau}_{o}^{w}& 0\end{array}\right|$

Finally,

$\sigma =\left|\begin{array}{cc}0& {\tau}_{i}^{w}\\ {\tau}_{o}^{w}& 0\end{array}\right|$ (7)

This is the deviatoric stress tensor, which describes the maximum shear forces experienced by an element in the core of the fluid. When the fluid is immovable this tensor is zero everywhere in the fluid. Hence, it is directly related to the velocity of the fluid.

It is well known that there are some tensor invariants like the trace and the determinant. Herein, a definition is proposed,

$\frac{1}{2}\rho {V}^{2}=\frac{1}{f}\sqrt{\left|\Vert \sigma \Vert \right|}$ (8)

*i*.*e*., the kinetic energy is proportional to the square root of the absolute value of the tensor’s
$\sigma $ determinant
$\Vert \sigma \Vert $. Practically, Equation (8) connects two invariant quantities for a given pressure gradient and geometry for a non-compressible fluid.

The coefficient 1/*f* identifies to the Fanning’s (not Darcy’s) friction factor and guarantees the steady state. Boundary conditions are satisfied given that at zero velocity the right-hand side of Equation (8) is zero, as shear stresses appear only and only if the fluid is in motion.

Before proceeding with the annuli problem, which illustrates the present methodology, consider a simple case, laminar flow in a cylindrical tube. Start by rearranging Equation (8),

${V}^{2}=\frac{2}{\rho f}\sqrt{\left|\Vert \sigma \Vert \right|}$ (9)

For this case,*
$\sigma $ * is symmetric. Thus,

$\sigma =\left|\begin{array}{cc}0& {\tau}^{w}\\ {\tau}^{w}& 0\end{array}\right|$ (10)

calculating the determinant,

$\sqrt{\left|\Vert \sigma \Vert \right|}=\sqrt{\left|0\times 0-{\tau}^{w}\times {\tau}^{w}\right|}={\tau}^{w}=\frac{\nabla p}{2}R\text{\hspace{0.17em}}\left(\text{N}/{\text{m}}^{\text{2}}\right)$ (11)

*Same time*, for laminar flow, from the empirical equation,

$\frac{1}{f}=\frac{Re}{16}$ (12)

Combining Equation (9), Equation (11) and Equation (12),

${V}^{2}=\frac{2}{\rho f}\frac{\nabla p}{2}R=\frac{2Re}{\rho 16}\frac{\nabla p}{2}R=\frac{2\rho VD}{\rho 16\mu}\frac{\nabla p}{2}R$ $\Rightarrow $

$V=\frac{D}{8\mu}\frac{\nabla p}{2}R$ (13)

Attention is needed to the fact that diameter *D* in Equation (12) is defined through the *Re* number and has a precise meaning. It is the quantity, which uniquely defines the area of the flow. In this case, the circular area of the tube is well defined by *D* (in order to define, for instance through *R*, we would need more parameters like 2*R* etc. For other less simple shapes it is trickier to define *D*).

Finally,

$V=\frac{\nabla p}{8\mu}{R}^{2}$ (14)

This last equation identifies with the well-known, analytically obtained, equation for the average velocity in a circular tube in laminar flow.

In a similar way the following well known equation was obtained which will be referred hereafter as the “classic one” for laminar flow,

$\frac{1}{2}\rho {V}^{2}=\frac{1}{f}4{\stackrel{\u02dc}{\tau}}_{w}$ with*
$\frac{1}{f}=\frac{Re}{64}$ *

${\stackrel{\u02dc}{\tau}}_{w}\equiv {\tau}_{\text{wallaverage}}$

and defined by the following equation,

${\stackrel{\u02dc}{\tau}}_{w}\left(2\pi R+2\pi KR\right)L=\Delta p\pi \left({R}^{2}-{K}^{2}{R}^{2}\right)$

${\stackrel{\u02dc}{\tau}}_{w}=\nabla p\frac{\pi \left({R}^{2}-{K}^{2}{R}^{2}\right)}{2\pi R+2\pi KR}=\nabla p\frac{A}{{P}_{\text{wetted}}}$

$4{\stackrel{\u02dc}{\tau}}_{w}=\nabla p\frac{4A}{{P}_{\text{wetted}}}=\nabla p{D}_{h}$

Hence,

$V=\frac{2{D}_{h}^{2}}{64\mu}\nabla p$ (15)

To find the diagonal matrix of the symmetric tensor Equation (10).

Let, $R=0.05\text{\hspace{0.17em}}\text{m}$, $\nabla p=0.1\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$, From (6) and $K=0$,

${\tau}^{w}=\frac{\nabla p}{2}R=0.0025\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

$\sigma =\left|\begin{array}{cc}0& {\tau}^{w}\\ {\tau}^{w}& 0\end{array}\right|=\left|\begin{array}{cc}0& 0.0025\\ 0.0025& 0\end{array}\right|$ (16)

and ${\sigma}_{\text{diagonal}}=\left|\begin{array}{cc}{\lambda}_{1}& 0\\ 0& {\lambda}_{2}\end{array}\right|$

where ${\lambda}_{1}$ and ${\lambda}_{2}$ are the two eigenvalues of $\sigma $ and thus,

${\sigma}_{\text{diagonal}}=\left|\begin{array}{cc}0.0025& 0\\ 0& -0.0025\end{array}\right|$ (17)

and the eigenvectors are,

*
${v}_{1}=\left|\text{\hspace{0.05em}}\begin{array}{c}1\\ 1\end{array}\text{\hspace{0.05em}}\right|$ * and
${v}_{2}=\left|\begin{array}{c}-1\\ 1\end{array}\right|$ (18)

${\mathrm{tan}}^{-1}1/4=45\u02da$

so, $\theta =45\u02da$.

In this case what happened is a rotation of an element under stresses that it is described by the stress tensor in Equation (16) and by an angle $\theta =45\u02da$. In this new virtual position, the element suffers no more shear stresses, just tensile stresses which are normal to the surface of the element.

Clearly shear stresses are due to the relative velocity, which exists in the laminar flow between adjacent elements, somehow proportional to the velocity as described in Equation (8).

*However now the eigenvector indicates a new direction*,*at *45˚*where shear stresses should disappear*.*The plains where shear stresses disappear are called *“*Principal-Plains*”.*What does it mean*?

It is clear now from Figure 2 what is happening. The eigenvector tells an element at the point *A*, “If you run faster than your neighbor bearing 45˚ to the left than you suffer no more shear stresses and you arrive to the point *C* where shear stresses do not exist as well”.

At point *C* the two neighboring elements are running with almost the same velocity (stationary point), thus shear stresses disappear. The situation is the same at point *A*, where the velocity is zero. Thus, line *AC* unites all points with zero shear stresses.

What is the gain finding a zero-shear stress direction?

Back to the less trivial situation, which reigns in an annular flow with $K=0.3$.

The parameters are (see Figure 3)

$R=0.05\text{\hspace{0.17em}}\text{m}$,*
$KR=0.015\text{\hspace{0.17em}}\text{m}$ *,

$\nabla p=0.1\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$,
$\lambda =0.615$,*
$\lambda R=0.03075\text{\hspace{0.17em}}\text{m}$ *

Calculating,

${\tau}_{i}^{w}=0.0024\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$ and ${\tau}_{o}^{w}=0.00155\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

with the help of Equation (5) and Equation (6),

$\sigma =\left|\begin{array}{cc}0& {\tau}_{i}^{w}\\ {\tau}_{o}^{w}& 0\end{array}\right|=\left|\begin{array}{cc}0& 0.0024\\ 0.00155& 0\end{array}\right|$ (19)

Figure 2. The velocity distribution for laminar flow in a cylindrical cross-sectional pipe.

Figure 3. Annular flow $K=0.3$, $R=0.05\text{\hspace{0.17em}}\text{m}$.

The eigenvalues are,*
${\lambda}_{1}=0.00193\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$ * and
${\lambda}_{2}=-0.00193\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

${\sigma}_{\text{diagonal}}=\left|\begin{array}{cc}0.00193& 0\\ 0& -0.00193\end{array}\right|$ (20)

and the eigenvectors are,

*
${v}_{1}=\left|\begin{array}{c}1.244\\ 1\end{array}\right|$ *and *
${v}_{2}=\left|\begin{array}{c}-1.244\\ 1\end{array}\right|$ * (21)

Leading to an angle ${\mathrm{tan}}^{-1}1.244/1=51.2\u02da$. So, $\theta =51.2\u02da$.

*Although the tensor in *Equation (19)*is not symmetric*,*the eigenvectors do exist **but are not orthogonal*.*In fact*,*they form an angle of
$2\theta =2\times {51.2}^{\circ}={102.4}^{\circ}$ *(see also Appendix B).

From Figure 3, it is clear that in this case the eigenvectors indicate a different direction than in the symmetric one, which is precisely the point *F* at a distance of 0.03075 m from the central axes Z.

In fact, there is a geometric relation between the angles and the parameter $\lambda $.

Following a rough approach, the point *A* sees
$\stackrel{\xaf}{CB}$ at an angle of 45˚ and respectively
$\stackrel{\xaf}{BF}$ at an angle of 51.2˚

$\frac{{51.2}^{\circ}}{{45}^{\circ}}=\frac{\stackrel{\xaf}{BF}}{\stackrel{\xaf}{CB}}=\frac{\stackrel{\xaf}{BF}}{0.0175}$

Hence, $\stackrel{\xaf}{BF}=0.0199$

$\stackrel{\xaf}{OF}=0.05-0.0199=0.0301$

However, $\stackrel{\xaf}{OF}=\lambda R$. Therefore, $\lambda \approx 0.6$.

In other words, the eigenvector indicates a direction to a point with cardinal significance for this kind of problems *i*.*e*., the point at maximum velocity.

*It is worth pointing out that as
$K\to 1$ *,*this very rough geometric approach used to identify parameter
$\lambda $ *, *becomes accurate*.*Instead*,*for
$K\to 0.0$ *,*there is lack of accuracy most likely due to the increasing asymmetry of the stresses*.*In Appendix B*,*there is a precise geometrical approach that demonstrates the relation of the angle*
$\theta $ *and the point of maximum velocity*.* Basically*,* it turns out that*
$\theta $ * is the average of the complementary angles the two eigenvectors produce*.

At this point, it is useful to produce an approximative geometric equation that relates the angle
$\theta $ with the parameters*
$\lambda $ * and*K*.

From Figure 3,

$\frac{\Theta}{{45}^{\circ}}=\frac{\stackrel{\xaf}{BF}}{\stackrel{\xaf}{CB}}$

and

*
$\stackrel{\xaf}{CB}=\frac{R-KR}{2}=\frac{R\left(1-K\right)}{2}$ *,
$\stackrel{\xaf}{BF}=R\left(1-\lambda \right)$

Finally,

$\lambda =1-\frac{\Theta}{{45}^{\circ}}\frac{1-K}{2}$ (22)

Equation (22) calculates $\lambda $ with good accuracy particularly for $0.3\le K\le 0.9$.

Moreover, Equation (22) demonstrates that $\lambda $ is somehow related to $\theta $, a property which will be useful later in the turbulent flow.

At this point, to calculate the average velocity for this annulus, the proposed Equation (8) for a supposed laminar flow is used,

Given,

$\frac{1}{2}\rho {V}^{2}=\frac{1}{f}\sqrt{\left|\Vert \sigma \Vert \right|}$ and $\frac{1}{f}=\frac{Re}{16}$

From the above two,

$V=\frac{D}{8\mu}\sqrt{\left|\Vert \sigma \Vert \right|}\text{\hspace{0.17em}}\left(\text{m}\cdot {\mathrm{sec}}^{-1}\right)$

An issue appears here; how to define*D*. Remember *D* comes from the definition of the *Re* and one must be very careful. A straightforward approach says that here the significant distance is the clearance of the annulus (the gap) *i*.*e*.,
$\left(R-KR\right)$. However, on the *r* axis and by symmetry (see Figure 1) there are two gaps.

So, ${D}_{h}=2\left(R-KR\right)=0.07\text{\hspace{0.17em}}\text{m}$, ${D}_{h}$ is called the Hydraulic Diameter.

On the other hand, the standard approach is to define the Hydraulic Diameter as,

${D}_{h}=\frac{4A}{P}\text{\hspace{0.17em}}(\; m\; )$

with*P* the wetted perimeter and *A* the cross-section area of the flow.

Fortunately, in this case above equation gives the same result, so ${D}_{h}=0.07\text{\hspace{0.17em}}\text{m}$ is accepted.

The invariant determinant of either matrix Equation (19) or Equation (20) is,

$\Vert \sigma \Vert =\left(-0.0024\times 0.00155\right)=\left(-0.00193\times 0.00193\right)$

Finally, $\sqrt{\left|\Vert \sigma \Vert \right|}={\lambda}_{1}\text{\hspace{0.17em}}\text{or}\text{\hspace{0.17em}}{\lambda}_{2}=0.00193\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$.

From Equation (21),

$V=\frac{{D}_{h}}{8\mu}\sqrt{\left|\Vert \sigma \Vert \right|}=\frac{0.07}{8\times {10}^{-3}}\times 0.00193=0.017\text{\hspace{0.17em}}\text{m}/\text{sec}$ (23)

From the exact analytical solution [7],

$V=\frac{\nabla p}{8\mu}{R}^{2}\left(\frac{1-{K}^{4}}{1-{K}^{2}}-\frac{1-{K}^{2}}{\mathrm{ln}\frac{1}{K}}\right)=\frac{0.1}{8\times {10}^{-3}}\times {0.05}^{2}\times 0.33=0.01\text{\hspace{0.17em}}\text{m}/\text{sec}$ (24)

Also, from the classical approach, Equation (15),

$V=\frac{2{D}_{h}^{2}}{64\mu}\nabla p=\frac{2\times {0.07}^{2}}{64\times {10}^{-3}}\times 0.1=0.0153\text{\hspace{0.17em}}\text{m}/\text{sec}$

Therefore, although the classical approach gives closer results to the analytical solution than Equation (23), they are still higher. There is need to investigate what is the issue here.

The first incorrect step was the use of the experimental Equation (12) coming from the measurements on cylindrical cross section pipes. The very definition of the Reynolds number demands geometric similarity. In fact, we cannot claim geometric similarity between an annulus and a circular pipe. The second incorrect step was the very fact that due to the non-symmetry of the wall stresses the velocity profile is skew. The simple use of the Hydraulic Diameter concept seems not to be sufficient to cure either problem. Practically this is exactly what the present work claims.

At this point, the definition of the traction vector responsible for the fluid’s motion is recalled,

$T=\left|\begin{array}{c}{T}_{r}\\ {T}_{z}\end{array}\right|=\sigma \cdot n=\lambda n$ (25)

In this case,
$n$ is the normalized eigenvector that defines the surface on which the traction vector applies. It is understood that on the principal plane the traction vector*
$T$ * and the eigenvector
$n$ are aligned.

The eigenvectors already provided this vector, however not in a normalized form. Therefore, to normalize it just remember that the angle $\theta $ must be 51.2˚.

Moreover, there is the trigonometric identity, $\mathrm{cos}{\theta}^{2}+\mathrm{sin}{\theta}^{2}=1$ so, inserting $\theta ={51.2}^{\circ}$

${0.627}^{2}+{0.779}^{2}=1$ (26)

$T=\left|\begin{array}{c}{T}_{r}\\ {T}_{z}\end{array}\right|=\left|\begin{array}{cc}0& 0.0024\\ 0.00155& 0\end{array}\right|\times \left|\begin{array}{c}0.627\\ 0.779\end{array}\right|=0.00193\times \left|\begin{array}{c}0.627\\ 0.779\end{array}\right|$

$\left|\begin{array}{c}{T}_{r}\\ {T}_{z}\end{array}\right|=\left|\begin{array}{c}0.00193\times 0.627\\ 0.00193\times 0.779\end{array}\right|=\left|\begin{array}{c}0.00121\\ 0.0015\end{array}\right|$

Given that $\theta ={51.2}^{\circ}>{45}^{\circ}$, the smallest component should identify with the ${T}_{z}$.

*To remember that the components of the traction vector are*,* by definition*,* orthogonal*,* and parallel to the relative axes*.

Finally,

$\left|\begin{array}{c}{T}_{z}\\ {T}_{r}\end{array}\right|=\left|\begin{array}{c}0.00193\times 0.627\\ 0.00193\times 0.779\end{array}\right|=\left|\begin{array}{c}0.00121\\ 0.0015\end{array}\right|$

Coming back to Equation (22),

$V=\frac{{D}_{h}}{8\mu}\sqrt{\left|\Vert \sigma \Vert \right|}=\frac{0.07}{8\times {10}^{-3}}\times 0.00193\times 0.627=0.0106\text{\hspace{0.17em}}\text{m}/\text{sec}$

An excellent fit with the analytical one.

Therefore, the equation proposed for a laminar flow should be,

$V=\frac{{D}_{h}}{8\mu}\sqrt{\left|\Vert \sigma \Vert \right|}\mathrm{cos}\theta $ (27)

Or in a more explicit way,

$V=\frac{{D}_{h}}{8\mu}\lambda \mathrm{cos}\theta \text{\hspace{0.05em}}\nabla p$ (28)

*I*.*e*., just the component along the z axis of the traction vector maters.

Equation (27)*gives an excellent correlation with the analytical *Equation (24),*for the whole range of the parameter K as we can observe from *Table 1 (A more extensive theoretical explanation is given in Appendix B).*It is now evident that the Hydraulic Diameter concept cures satisfactorily the similarity problem*,*however*, *the main problem due to the skew profile is cured by considering just the Z component of the traction vector*.

In Table 1, it is evident that Equation (15) fails consistently. This is reasonable due to the fact that it is based on the semiempirical Equation (12), which in turn

Table 1. Average velocities m/sec for different *K* parameters and for laminar flow.

is based on values obtained in circular pipes. This problem is supposedly cured given the adoption of the Hydraulic Diameter, instead it is not. We believe that the cause of the insisting discrepancies is the skew velocity profile due to the non-symmetry of the flow.

For this reason, the correction factor ${\Phi}^{*}$ to Equation (15) as suggested and defined in [4] is introduced (see Appendix B),

$\frac{{D}_{eff}}{{D}_{h}}={\Phi}^{*}$ (29)

*It seems that in * *[4]* (*Jones et al.*)*practically tried to correct the problem arising from the lack of geometric similarity and the skewness*,*as explained*,*modifying Reynolds number by reducing the Hydraulic Diameter of the flow*.*Instead*,*in the present method just the component of the traction force along the z axis has been considered*.*Finally*,*it appears that this very component is responsible for the motion of the fluid forwards*.*Hence*,*two distinct approaches are recognized here*.

Moreover, in Table 1, it has been included the average velocity, calculated as “Flow in a duct” by the well-known and analytically obtained equation,

$V=\frac{{h}^{2}}{3\mu}\nabla p$ (30)

Here *h* equals half the clearance *i*.*e*.,

$h=\frac{1}{2}\left(R-KR\right)$

Equation (30)* is obtained on the assumption that the clearance *2*h is substantially smaller than the wide of the duct*.* It is interesting to observe that for
$K\ge 0.3$ *, Equation (30)* provides the most reliable and easy way to calculate the average velocity for flow through concentric annuli*,* at least for laminar flow*.* Therefore*,* a promising idea could be to transform a *“*duct problem*”* into an annular flow problem with
$K\ge 0.8$ and vice-versa*.*To see how the above *Equation (27)*works for laminar flow *(see Appendix E).

2.2. Flow through an Annulus-Turbulent Flow

To continue with the illustrative examples to better understand this new method, consider the case of turbulent flow. As it is mentioned in laminar flow, the concept of average velocity is a mathematical one. On the other hand, in the developed turbulent flow this concept has, along with the mathematical meaning, a physical one. The reason for this is that in the so-called “plug flow” the velocity is quasi constant at any point in the core of the fluid, except of course very close to the walls.

So, in a similar way with Equation (8),

$\frac{1}{2}\rho {V}^{2}=\frac{1}{f}\sqrt{\left|\Vert \sigma \Vert \right|}$ (31)

with the parameter *f* to be function of,

$f=F\left(V,{D}_{h},\mu ,\rho ,\epsilon \right)=F\left(R{e}_{h},\epsilon \right)$ (32)

At the same time, there is the classical definition,

$\frac{1}{2}\rho {V}^{2}=\frac{4{\tau}_{w}}{f}=\frac{4}{f}\frac{\nabla p}{4}D$ (33)

*In *Equation (33)*f* *stands for Darcy*’*s friction factor and it is *4*times greater than the f in *Equation (31)*that stands for Fanning*’*s factor*.

From Equation (33),

*
$V=\frac{1}{\sqrt{f}}\sqrt{\frac{2D\nabla p}{\rho}}$ * (34)

Now, in order to compare Equation (31) and Equation (34) in turbulent flow, the velocity should be fixed so to define the Reynolds number, *i*.*e*., we will solve the above equations for
$\nabla p$

From Equation (5) and Equation (6) we have,

${\tau}_{o}^{w}=\frac{\nabla p}{2}R\left(1-{\lambda}^{2}\right)$ (35)

${\tau}_{i}^{w}=\frac{\nabla p}{2}R\left(K-\frac{{\lambda}^{2}}{K}\right)$ (36)

In turbulent flow, the shear stresses in the core of the fluid due to cohesiveness do not exist. Under high pressure to achieve higher velocities the fluid lost cohesiveness. It is not any more considered a continuum according to the meaning we give in mechanics. Shear stresses, such as these defined in laminar flow, due to the relative velocities between adjacent fluid elements, do not exist. The fluid is in a new and different state. Looking on Moody’s diagram, the sharp discontinuity between laminar and turbulent flow is evident. Imagine the transformation happening to the soil under high stresses during a strong earthquake. The soil fluidifies and behaves like moving sand. The fluid in turbulent flow behaves like a cluster of moving particles (the eddies) with high velocity and friction occurring from the collisions with each other and the wall. A similar situation exists during heat transfer in turbulent flow. In this case, the heat transfer is almost exclusively due to convection and not to conduction, which needs the concept of a mutual contact.

Finally, the stress tensor becomes,

$\sigma =\left|\begin{array}{cc}0& {\tau}_{i}^{w}\\ {\tau}_{o}^{w}& 0\end{array}\right|=\frac{\nabla p}{2}R\left|\begin{array}{cc}0& K-\frac{{\lambda}^{2}}{K}\\ 1-{\lambda}^{2}& 0\end{array}\right|$ (37)

From Equation (37) it is more than clear that the orientation of the eigenvectors is now just a matter of geometry.

And the determinant,

$\Vert \sigma \Vert =-\frac{\nabla p}{2}R\left(1-{\lambda}^{2}\right)\frac{\nabla p}{2}R\left(K-\frac{{\lambda}^{2}}{K}\right)$

$\sqrt{\left|\Vert \sigma \Vert \right|}=\sqrt{\left|-\frac{\nabla p}{2}R\left(1-{\lambda}^{2}\right)\frac{\nabla p}{2}R\left(K-\frac{{\lambda}^{2}}{K}\right)\right|}=\frac{\nabla p}{2}R\sqrt{\left|\left(1-{\lambda}^{2}\right)\left(K-\frac{{\lambda}^{2}}{K}\right)\right|}$

Inserting this into Equation (31),

$\frac{1}{2}\rho {V}^{2}=\frac{1}{f}\frac{\nabla p}{2}R\sqrt{\left|\left(1-{\lambda}^{2}\right)\left(K-\frac{{\lambda}^{2}}{K}\right)\right|}$

Finally, Equation (31) becomes,

$V=2\frac{1}{\sqrt{f}}\sqrt{\frac{R\sqrt{\left|\left(1-{\lambda}^{2}\right)\left(K-\frac{{\lambda}^{2}}{K}\right)\right|}}{\rho}}\sqrt{\nabla p}$ (38)

*In *Equation (38),*the Darcy*’*s friction factor is now accepted*.*It can be obtained through *Equation (40)*or *Equation (1).*The extra factor *2*in *Equation (38)*takes care of the fact that the Fanning*’*s friction factor is *4*times smaller the Darcy*’*s one*.

Similarly, Equation (34) becomes,

$V=\frac{1}{\sqrt{f}}\sqrt{\frac{2D}{\rho}}\sqrt{\nabla p}$ (39)

Ready now to compare the two approaches for the case of $K=0.2$ annular flow of water with $Re={10}^{5}$.

Parameters given, $\rho =1000\text{\hspace{0.17em}}\text{kg}/{\text{m}}^{\text{3}}$

$\mu ={10}^{-3}\text{kg}/\text{m}\cdot \text{sec}$,*
$R=0.05\text{\hspace{0.17em}}\text{m}$ *

$\lambda =0.546$,*
${d}_{h}=0.08\text{\hspace{0.17em}}\text{m}$ *

$R{e}_{d}=\frac{\rho V{d}_{h}}{\mu}=\frac{1000\times V\times 0.08}{{10}^{-3}}={10}^{5}$ *
$\to $ *
$V=1.25\text{\hspace{0.17em}}\text{m}/\text{sec}$

with the help of Equation (35), Equation (36),

${\tau}_{i}^{w}=1.29\times R\frac{\nabla p}{2}$ and*
${\tau}_{o}^{w}=0.702\times R\frac{\nabla p}{2}$ *

Hence,

$\sigma =\left|\begin{array}{cc}0& {\tau}_{i}^{w}\\ {\tau}_{o}^{w}& 0\end{array}\right|=R\frac{\nabla p}{2}\left|\begin{array}{cc}0& 1.29\\ 0.702& 0\end{array}\right|$ $\to $

${\lambda}_{1}$ or*
${\lambda}_{2}=\pm 0.952\times R\frac{\nabla p}{2}$ *

*
${v}_{1}=\left|\begin{array}{c}1.356\\ 1\end{array}\right|$ * and*
${v}_{2}=\left|\begin{array}{c}-1.356\\ 1\end{array}\right|$ *

and, $\theta ={53.6}^{\circ}$.

To calculate the friction factor, the approximative but explicit, equation from Haaland [8] and for smooth pipes, is used,

$\frac{1}{\sqrt{f}}=-1.8\mathrm{log}\left(\frac{6.9}{R{e}_{d}}\right)$ (40)

Valid for, $4\times {10}^{4}<Re<{10}^{8}$

$\frac{1}{\sqrt{f}}=-1.8\mathrm{log}\left(\frac{6.9}{{10}^{5}}\right)=7.49$ (41)

Now ready to calculate the
$\nabla p$, firstly from Equation (38), then with the method proposed by Jones *et al.* [4], and finally compare.

$1.25=2\times 7.49\times \sqrt{\frac{2\times \sqrt{\left|\left(0.952\times R\frac{\nabla p}{2}\right)\left(-0.952\times R\frac{\nabla p}{2}\right)\right|}}{1000}}$

$1.25=2\times 7.49\times \sqrt{\frac{R\times 0.952}{1000}}\times \sqrt{\nabla p}$

$1.25=2\times 7.49\times \sqrt{\frac{0.05\times 0.952}{1000}}\times \sqrt{\nabla p}$
$\to $ *
$\nabla p=146\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$ * (42)

At this point following the route proposed (Jones *et al.*) [4], which has been generally accepted [2] (see Appendix B)

$R{e}_{d}^{\ast}=R{e}_{d}\times {\Phi}^{*}$

For
$K=0.2$ *
$\to $ *
${\Phi}^{*}=0.692$. Thus,*
$R{e}_{d}^{\ast}=69200$ *

$\frac{1}{\sqrt{f}}=2\mathrm{log}\left(R{e}_{d}^{\ast}\times \sqrt{f}\right)-0.8$ (43)

Guess,*
$f=0.0193$
$\to $ *
$7.2=7.916-0.8=7.17$

$\to $ *
$f=0.0195$ *

Enter,*
$f=0.0195$ **
$\to $ *
$7.16=7.97-0.8=7.17$

$\to $ *
$f=0.0195$ * O.K. accepted.

And from Equation (39),

$1.25=7.16\sqrt{\frac{2\times 0.08}{1000}}\times \sqrt{\nabla p}$ $\to $ $\nabla p=191\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$ (44)

So, there is here an unacceptable difference of almost 40%.

In Table 2, one can observe interesting outcomes.

- Equation (38) always gives lower results than Equation (39)*.

- The differences decrease as *Κ* parameter increases.

- Most probably the discrepancies are related to the geometry of the annulus. The discrepancies due to the non-geometric similarity are cured, given that the concept of Hydraulic Diameter has been used. Hence, the insisting discrepancies are to be attributed to the non-symmetry of the flow.

- As $K\to 0$, the annulus practically tends to become a standard cylindrical cross-section tube, not geometrically similar to an annulus. Hence, it is reasonable that Equation (38) is not behaving correctly.

- On the other hand, as $K\to 1$, the annulus tends to become a perfect duct, geometric non-similarity remains, however, skewness is reduced. This fact explains the convergence of Equation (38) to Equation (39)*.

Table 2. Pressure gradient in N/m^{3}.

- So, there is need to correct Equation (38), however, not in the Jones’ way; just by changing the Hydraulic Diameter, and then altering Reynolds number. Instead, correct by using a correction factor, which helps to reduce the problem due to the skewness of the velocity’s profile coming from the non-symmetric wall stresses.

- The correction factor should be strongly related to this new geometry, *i*.*e*., it should best describe the geometry.

- The correction factor should converge to zero as $K\to 0$

- The correction factor should converge to unity as
$K\to 1$, *i*.*e*., to the duct flow where skewness disappears.

So, the best candidate according to the Equation (22) seems to be $\lambda $. Moreover, it is this parameter that describes the peculiar geometry of the system to its fullest. In fact, from the approximative Equation (22),

$\lambda =1-\frac{\Theta}{{45}^{\circ}}\frac{1-K}{2}$ (22)

as,

$K\to 1$, $\lambda \to 1$ and $\Theta \to {45}^{\circ}$

$K\to 0$, $\lambda \to 0$ and $\Theta \to {90}^{\circ}$

So, Equation (31) could become,

$\frac{1}{2}\rho {V}^{2}=\frac{1}{f}\sqrt{\lambda \left|\Vert \sigma \Vert \right|}$

And, finally, from Equation (38),

$V=2\frac{1}{\sqrt{f}}\sqrt{\frac{R\sqrt{\left(1-{\lambda}^{2}\right)\left(K-\frac{{\lambda}^{2}}{K}\right)}}{\rho}}\sqrt{\sqrt{\lambda}}\times \sqrt{\nabla p}$ (45)

*For
$K\to 1$ and
$\lambda \to 1$ *,*we get
$V\to 0$ *,*which is acceptable*,*because the flow turns into duct flow and the duct*’*s gap is going to close*.*On the contrary*,*as
$K\to 0$ * Equation (45)*blows up*,*it is not any more well defined*.*This is a reasonable and a well expected outcome for *Equation (45),*i*.*e*., Equation (45)*doesn*’*t**work for
$K=0$ *.*It doesn*’*t make sense*.*The geometry is a perfect circular cross section pipe*,*and the standard equations are valid *(see Appendix D).*The matrix for
$K=0$ is not diagonalizable*.*Moreover*,*at
$K=0$ * *even**
$\lambda $ **does not make sense as a logarithmic average*.

So, in this case,*
$\sqrt{\sqrt{\lambda}}=\sqrt{\sqrt{0.546}}=0.86$ *

$1.25=2\times 7.49\times 0.0069\times 0.86\times \sqrt{\nabla p}$ $\to $ $\nabla p=198\left(\text{N}/{\text{m}}^{\text{3}}\right)$

From Table 3 interesting facts appear (see also Figure 4 and Figure 5).

$\nabla p$ -Equation (39) the classic and $\nabla p$ -Equation (39)* the classic* diverge as

Figure 4.
$\nabla p$ in N/m^{3} for (
$0.05\le K\le 0.6$ ). Data from Table 3.

Figure 5.
$\nabla p$ in N/m^{3} for (
$0.05\le K\le 0.8$ ). Data from Table 3.

Table 3. The pressure gradient
$\nabla p$ in N/m^{3}, turbulent flow in annuli (
$R{e}_{d}={10}^{5}$ ) and for
$0.05\le K\le 0.9$.

$K\to 0.9$, compelling the fact that the Prandl’s Equation (1), corrected with the help of the correction factor
${\Phi}^{*}$, [4] does not seems to hold good for ducts, *i*.*e*., where the non-symmetry is not the main issue.

Instead,
$\nabla p$ -Equation (39) the classic and
$\nabla p$ -Equation (38)
$\sqrt{\sqrt{\lambda}}$ converge quickly as
$K\to 0.9$ compelling the fact that the two approaches (*i*.*e*., Equation (39)* the classic* and Equation (38)
$\sqrt{\sqrt{\lambda}}$ ) are not equivalent.

$\nabla p$ -Equation (38)
$\sqrt{\sqrt{\lambda}}$ correlates very well with the generally acceptable data that the
$\nabla p$ -Equation (39)* *the classic* *provides, particularly for
$0.1\le K\le 0.9$. However,
$\nabla p$ -Equation (39)* *the classic** diverges slightly as
$K\to 0.9$.

*An interesting*,*yet weak correlation appears also with
$\nabla p$ -*Equation (38)
$\sqrt{\sqrt{\mathrm{cos}\Theta}}$.*In fact*,*just remembering that the tensor
$\sigma $ is built on the two eigenvectors*,*
$\mathrm{cos}\Theta $ **has to do with the dot product of these eigenvectors*.*In fact*,*the angle of these two vectors is exactly
$2\Theta $ *.*This fact explains the presence of the correlation coefficient as
$\sqrt{\sqrt{\lambda}}$ * (see also Appendix C).

3. About Rectangular Ducts

Obviously, rectangular ducts are not geometrically similar to a cylindrical cross section and the use of the Hydraulic Diameter concept is necessary. Moreover, it will be shown immediately that in this case it is also sufficient given the flow is symmetric. So, will be evident how the simple use of Hydraulic Diameter is sufficient to describe this kind of flow. These findings are in contrast with [3] (Jones *et al.*) and [2] that claim there is need for a correction factor function to the apex ratio of the duct.

3.1. Turbulent Flow in Rectangular Ducts

The findings about the concentric annuli, Equation (45) in particular, could be a useful platform to calculate friction factors for turbulent flow in any kind of rectangular ducts. For a laminar flow, there is a simple and precise analytical solution, however, this subject is coming later.

Considering the simplest case of a duct with a width much larger than the spacing, *i*.*e*.,
$w\gg s$.

Same time imagine an annulus with $K=0.9$, $R=10\text{\hspace{0.17em}}\text{m}$. In this case, the gap (spacing) $s=1\text{\hspace{0.17em}}\text{m}$ and $w=63\text{\hspace{0.17em}}\text{m}$. Therefore, the condition $w\gg s$ is almost met.

It is possible to compare the two approaches (the present method against [3] (Jones *et al.*), by calculating the pressure drop in turbulent flow with these two distinct methods. Suppose the parameters are,

$K=0.9$ annular flow of water
$V=0.5\text{\hspace{0.17em}}\text{m}/\text{sec}$,*
$\rho =1000\text{\hspace{0.17em}}\text{kg}/{\text{m}}^{\text{3}}$ *

$\mu ={10}^{-3}\text{kg}/\text{m}\cdot \mathrm{sec}$, $\lambda =0.95$, $R=10\text{\hspace{0.17em}}\text{m}$

$R{e}_{d}=\frac{\rho V{d}_{h}}{\mu}=\frac{1000\times 0.5\times 2}{{10}^{-3}}={10}^{6}$

From Equation (40),

$\frac{1}{\sqrt{f}}=-1.8\mathrm{log}\frac{6.9}{{10}^{6}}$
$\to $ *
$\frac{1}{\sqrt{f}}=9.29$ *

From Equation (45),

$0.5=2\times 9.29\sqrt{\frac{10\sqrt{\left(1-{0.95}^{2}\right)\left(0.9-\frac{{0.95}^{2}}{0.9}\right)}}{\rho}}\sqrt{\sqrt{0.95}}\times \sqrt{\nabla p}$ $\to $ $\nabla p=0.74\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$

Hence, using Equation (45) as a platform, a pressure gradient 0.74 N/m^{3} has been obtained for a duct with
$w\gg s$ and
$s=1$.

And now using the approach by [3] (Jones *et al.*) (see Appendix B, Equation (B8)).

For
$\frac{s}{w}=\frac{1}{63}$ *
$\to $ *
${\Phi}^{*}=0.68$

${d}_{h}=\frac{4\times w\times 1}{2w}=2\text{\hspace{0.17em}}\text{m}$

$R{e}_{{d}_{h}}\times {\Phi}^{*}=\frac{\rho V{d}_{h}}{\mu}{\Phi}^{*}=\frac{1000\times 0.5\times 2}{{10}^{-3}}\times 0.68$

$R{e}_{d}^{\ast}=680435$

$\frac{1}{\sqrt{f}}=-1.8\mathrm{log}\frac{6.9}{680435}=9$

From Equation (39),

$0.5=9\sqrt{\frac{2\times 2}{1000}}\times \sqrt{\nabla p}$

$\nabla p=0.77\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$

So, we observe an excellent fit.

Moreover, checking this approach for a real rectangular duct with $s/w=1/3$.

From Figure 6, it is clear what is the idea. The fact is that Equation (45) for concentric annuli with*
$K=0.9$ * could be a promising platform to find an exact solution for any rectangular duct in turbulent flow.

Obviously, any point in the annulus shares the same $\nabla p$ and this is also valid inside the element with $w=3\text{\hspace{0.17em}}\text{m}$ in Figure 6.

The similarity condition for the platform is met, given that the decision was to transform any rectangular duct to an annulus always with $K=0.9$. What is left is to adapt the geometry of the annulus to share the same ${d}_{h}$ with the duct.

Find the Hydraulic diameter of the duct,

${d}_{h}=\frac{4\times 3\times 1}{2\times 3+2\times 1}=\frac{12}{8}=1.5\text{\hspace{0.17em}}\text{m}$

For the annulus, it should be,

$2R\left(1-K\right)=1.5\text{\hspace{0.17em}}\text{m}$

$R=7.5\text{\hspace{0.17em}}\text{m}$

Figure 6. Annulus with $K=0.9$ as a platform to calculate turbulent flow in rectangular ducts.

*Remember*,*for the annuli*,*
${d}_{h}=2\times gap$ *

$R{e}_{{d}_{h}}=\frac{\rho V{d}_{h}}{\mu}=\frac{1000\times 0.5\times 1.5}{{10}^{-3}}=7.5\times {10}^{5}$

*And from *Equation (40)
$\frac{1}{\sqrt{f}}=9.06$.

Hence, from Equation (45),

$0.5=2\times 9.06\sqrt{\frac{7.5\sqrt{\left(1-{0.95}^{2}\right)\left(0.9-\frac{{0.95}^{2}}{0.9}\right)}}{1000}}\sqrt{\sqrt{0.95}}\times \sqrt{\nabla p}$

$\nabla p=1.04\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$

*Similarly*,*in* [3] (*Jones et al.*) (see Appendix B, Equation (B8)).

$R{e}_{{d}_{h}}^{\ast}=\frac{\rho V{d}_{h}}{\mu}{\Phi}^{*}=\frac{1000\times 0.5\times 1.5}{{10}^{-3}}\times 0.92=6.9\times {10}^{5}$

$\frac{1}{\sqrt{f}}=9$
$\to $ *
$0.5=9\sqrt{\frac{2\times 1.5}{1000}}\times \sqrt{\nabla p}$ *
$\nabla p=1.03\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$

Again, an excellent fit.

Interesting facts appear in Table 4:

- Seems now possible to calculate with excellent accuracy the pressure drops in any rectangular duct using as a platform an annular geometry with $K=0.9$ and a radius, such as both Hydraulic Diameters be the same.

- The correlation between the two methods is excellent, except for the case of the square ducts where we see 4.5% difference. It is possible to improve this by eliminating the curvature factor $\sqrt{\sqrt{\lambda}}$ from Equation (45) in all the cases, which deal with square ducts.

- Moreover, the correlation with [4] (Jones *et al.*), for the case
$\frac{s}{w}\to 0$ or

$s\ll w$ is not good. This fact is strange, given it is under such a condition the duct tent to identify with an annulus. An applaudable explanation could be that in this case the flow is quasi symmetrical, and the method proposed in [3] is not behaving well.

- Notice, that although the two-aspect ratios 2/2 and 1/4 are different still have same cross section area. Nevertheless, they give different pressure drop. This is exactly the reason why we need the concept of the Hydraulic Diameter.

- Furthermore, the square geometry 2/2 gives lower pressure drop than 1/4, because the wetted perimeter is smaller. The circular cross section tube, which Colebrook used, will provide the lower pressure drop relative to the involved cross section area. Hence, it is evident that the circular cross section is the best choice to design an efficient hydraulic conduit.

Simplifying: having defined as platform an annulus with
$K=0.9$ with the help of the following equation, the *R *of the annulus is obtained.

$2R\left(1-K\right)={d}_{h}$

Table 4.
$\nabla p$ in N/m^{3} turbulent flow in rectangular ducts with different aspect ratio. Equation (45) is used with
$K=0.9$.

In this case, ${d}_{h}$ is the Hydraulic Diameter of the duct under consideration and since $K=0.9$ is always maintained the above equation becomes,

$R=\frac{{d}_{h}}{0.2}$

Thus, the Equation (45) suggested as platform to solve any duct problem, will simply be,

$V=2\frac{1}{\sqrt{f}}\sqrt{\frac{R\times 0.1}{\rho}}\times 0.987\times \sqrt{\nabla p}$

or

$V=\frac{1}{\sqrt{f}}\sqrt{\frac{2{d}_{h}}{\rho}}\times 0.987\times \sqrt{\nabla p}$ (46)

Valid for any rectangular duct in turbulent flow.

Practically, Equation (46) almost identifies with Equation (39) for cylindrical cross section pipes in turbulent flow, except for the coefficient 0.987.

This remarkable finding demonstrates that the use of the Hydraulic Diameter is sufficient to correlate Equation (39) to real rectangular ducts problems. Both equations, Equation (46) and Equation (39), depend only to the ${d}_{h}$ and nothing more. In fact, this is compelling evidence that the skewness of the profile matters more than the lack of geometric similarity. In this case, skewness is absent. Therefore, it seems that there is no need to recur to the coefficient ${\Phi}^{*}$ [2] [3] or further complications, at least for rectangular ducts.

And for a square duct, in turbulent flow,

$V=\frac{1}{\sqrt{f}}\sqrt{\frac{2{d}_{h}}{\rho}}\times \sqrt{\nabla p}$ (47)

with ${d}_{h}\equiv side$.

Solving by this method the example 6.15 on page 344 in [7]

It is a square duct 0.25 m × 0.25 m

The fluid is air*
$\rho =1.22\text{\hspace{0.17em}}\text{kg}/{\text{m}}^{\text{3}}$ *,
$\nu =1.45\times {10}^{-5}{\text{m}}^{\text{2}}/\text{sec}$

$L=30\text{\hspace{0.17em}}\text{m}$,
$V=11.2\text{\hspace{0.17em}}\text{m}/\text{sec}$,*
$\epsilon =1\times {10}^{-4}\text{\hspace{0.05em}}\text{\hspace{0.05em}}\text{m}$ *

Find the pressure drop*
$\Delta p$ *in N/m^{2}

First find the*
${d}_{h}$ *

For this square duct ${d}_{h}=0.25\text{\hspace{0.17em}}\text{m}$

Now, an annulus with the same*
${d}_{h}$ given *that the*
$K=0.9$ should be defined*,

Thus, $2R\left(1-K\right)=0.25\text{\hspace{0.17em}}\text{m}$ $\to $ $R=1.25\text{\hspace{0.17em}}\text{m}$

Also, the gap = 0.125 m

Finding Reynolds number,

$R{e}_{{d}_{h}}=\frac{11.2\times 0.25}{1.45\times {10}^{-5}}=193100$

Using the Haaland formula [8],

$\frac{1}{\sqrt{f}}=-1.8\mathrm{log}\left\{\frac{6.9}{193100}+{\left(\frac{{10}^{-4}/0.25}{3.7}\right)}^{1.11}\right\}$

$\frac{1}{\sqrt{f}}=7.42$

Hence, from Equation (47),

$11.2=7.42\sqrt{\frac{2\times 0.25}{1,22}}\sqrt{\nabla p}$

$\nabla p=5.5\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$ or $\Delta p=165\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

From the book, $\Delta p=163\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$.

3.2. Laminar Flow in Rectangular Ducts

For laminar flow in rectangular ducts exact analytical solutions are available. The general Equation (27) correlates quite well with the analytic solutions.

$V=\frac{{D}_{h}}{8\mu}\sqrt{\left|\Vert \sigma \Vert \right|}\mathrm{cos}\theta $ (27)

It is not claimed here that using Equation (27) is more practical than the known analytical equations, however, it helps to view the problem under a different perspective and understand better the dynamics of the flow. On the other hand Equation (27) is useful in cases the flow in a conduit is skewed for other reasons (see paragraph 4).

Therefore, Equation (27) for laminar flow in a duct becomes,

$V=\frac{{D}_{h}}{8\mu}{\tau}_{w}\mathrm{cos}\theta =\frac{{D}_{h}}{8\mu}h\nabla pcos\theta $

Given that, for the conjugate annulus, $K=0.9$ then $\mathrm{cos}\theta =0.698$

$V=\frac{{D}_{h}}{8\mu}\sqrt{\left|\Vert \sigma \Vert \right|}\times 0.698$ (27a)

$V=\frac{{D}_{h}}{8\mu}h\nabla p\times 0.698$ (27b)

Both equations are valid for any rectangular duct in laminar flow.

In this case, *h* is half the gap of the duct. It is claimed here that (27b) is more reliable than the exact analytic equation for a duct, *i*.*e*., Equation (30). Practically, although Equation (27b) almost identifies with Equation (30), the former has been obtained without the condition
$s\ll w$.

Indeed,

$V=\frac{{D}_{h}}{8\mu}h\nabla p\times 0.698=\frac{4h}{8\mu}h\nabla p\times 0.698$

$V=\frac{{h}^{2}}{2\mu}\nabla p\times 0.698~\frac{{h}^{2}}{3\mu}\nabla p\text{\hspace{0.17em}}\text{m}/\text{s}$

4. Other Findings

Another case of a non-symmetrical flow is the flow of two adjacent immiscible fluids in a duct [7] [9]. However, in this case, the non-symmetry comes from the differences on the viscosities, instead of the geometry. It would be interesting to observe how the present method could figure out the average velocities in this classic problem.

The expected velocity distribution at a steady state for the two fluids is represented by the curve $\left(\alpha ,\beta ,\gamma ,\delta \right)$.

At the point, $x=L$, the shear stresses, disappear, point of maximum velocity.

The first step is to calculate the stresses on the walls.

$\nabla \tau =\nabla p$ $\to $ $\text{d}\tau =\nabla p\text{d}x$ $\to $

${\tau}^{\prime}=\nabla p\times x+{C}^{\prime}$ and ${\tau}^{\u2033}=\nabla p\times x+{C}^{\u2033}$

At $x=0$ $\to $ ${\tau}^{\prime}={\tau}^{\u2033}$ and so,

${C}^{\prime}={C}^{\u2033}=-\nabla p\times L$

Finally,

${\tau}^{\prime}=\nabla p\left(x-L\right)$ and ${\tau}^{\u2033}=\nabla p\left(x-L\right)$

So,

${{\tau}^{\prime}}_{w}=-\nabla p\left(b+L\right)$

*
${{\tau}^{\u2033}}_{w}=\nabla p\left(b-L\right)$ * (47)

Now,

${{\tau}^{\prime}}_{w}-{{\tau}^{\u2033}}_{w}=-\nabla p\times 2b$ (48)

${{\tau}^{\prime}}_{w}+{{\tau}^{\u2033}}_{w}=\nabla p\times 2L$ (49)

Hence, there are two equations in two unknowns and a parameter*L*.

At the same time,

$\left|{{\tau}^{\prime}}_{w}\right|-\left|{{\tau}^{\u2033}}_{w}\right|=\nabla p\times 2L$ (50)

From Equation (50), it is evident that *L* is proportional to the difference of the absolute values of the wall stresses, *i*.*e*., *L* is proportional to the non-symmetry of the flow therefor, to the difference of the viscosities.

$L\propto \left({\mu}^{\prime}-{\mu}^{\u2033}\right)$

Or to get rid of the problem with the units,

$\frac{L}{b}\propto \frac{{\mu}^{\prime}-{\mu}^{\u2033}}{{\mu}^{\prime}+{\mu}^{\u2033}}$

However, for ${\mu}^{\prime}\gg {\mu}^{\u2033}$ the flow should become a normal laminar flow in a rectangular duct with $L=b/2$.

Finally,

$L=\frac{1}{2}\times \frac{{\mu}^{\prime}-{\mu}^{\u2033}}{{\mu}^{\prime}+{\mu}^{\u2033}}b\text{\hspace{0.17em}}\left(\text{m}\right)$ (51)

Let now,

${\mu}^{\prime}=1\times {10}^{-3}\text{kg}/\text{m}\cdot \mathrm{sec}$

${\mu}^{\u2033}=0.5\times {10}^{-3}\text{kg}/\text{m}\cdot \mathrm{sec}$

$b=0.025\text{\hspace{0.17em}}\text{m}$, $\nabla p=0.1\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$

From Equation (51), $L=0.0042\text{\hspace{0.17em}}\text{m}$.

For to find the average velocity
${V}^{\u2033}\left(\text{m}/\text{sec}\right)$ focus on the curve (*α*, *β*) in Figure 7. Imagine that the fluid with
${\mu}^{\u2033}$ is moving in a rectangular duct with a gap equal to *b*. So, Equation (27a) can be used even though the flow is skew.

$V=\frac{{D}_{h}}{8\mu}\sqrt{\left|\Vert \sigma \Vert \right|}\times 0.698$ (27a)

However, to proceed it is necessary to build the stress tensor,

$\sigma =\left|\begin{array}{cc}0& {{\tau}^{\u2033}}_{w}\\ {{\tau}^{\u2033}}_{?}& 0\end{array}\right|$

the stress on the wall it is known,

${{\tau}^{\u2033}}_{w}=\nabla p\left(b-L\right)=0.00208\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

However, there is need to know the conjugate stress, which provides the rest of the velocity distribution (*α*, *β*) in Figure 7. To figure out this stress, there are two ways. Either to imagine that the fluid has a skew distribution (*α*, *β*, *ζ*) or just to increase the space and imagine a symmetric laminar flow with new symmetric and an almost parabolic distribution (*α*,* β*, *ξ*) (see Figure 8).

First route,

$\mathrm{tan}\theta =\frac{0.025-0.0042}{0.0125}=1.664$

$\theta ={59}^{\circ}$

From Equation (E.2),

$\mathrm{tan}2\left(90-\theta \right)=\frac{\sigma}{\frac{{\tau}_{12}-{\tau}_{21}}{2}}$ (E.2)

$\mathrm{tan}2\left(90-59\right)=2\frac{\sigma}{{\tau}_{12}-{\tau}_{21}}=2\frac{\sqrt{{{\tau}^{\u2033}}_{w}\times {{\tau}^{\u2033}}_{?}}}{{\tau}_{12}-{\tau}_{21}}$

Or $1.88=2\frac{\sqrt{0.00208\times {{\tau}^{\u2033}}_{?}}}{0.00208-{{\tau}^{\u2033}}_{w}}$ obviously two values for ${{\tau}^{\u2033}}_{?}$ are provided,

${{\tau}^{\u2033}}_{?}=0.00576\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$ and*
${{\tau}^{\u2033}}_{?}=0.000746\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$ *

Figure 7. Flow of two immiscible fluids with different viscosities ${\mu}^{\prime},{\mu}^{\u2033}\left(\text{kg}/\text{m}\cdot \mathrm{sec}\right)$ with ${\mu}^{\prime}>{\mu}^{\u2033}$.

Figure 8. Flow of two immiscible fluids with different viscosities ${\mu}^{\prime},{\mu}^{\u2033}\left(\text{kg}/\text{m}\cdot \mathrm{sec}\right)$.

The correct value should be greater than 0.00208 N/m^{2},

$\sigma =\left|\begin{array}{cc}0& 0.00576\\ 0.00208& 0\end{array}\right|$

$\lambda =0.00346$ and ${v}_{1}=\left|\begin{array}{c}1.664\\ 1\end{array}\right|$ or $\theta ={59}^{\circ}$

From Equation (27a) with*
${D}_{h}=2\text{\hspace{0.17em}}times\text{\hspace{0.17em}}b$ *

${V}^{\u2033}=\frac{2\times 0.025}{8\times 0.5\times {10}^{-3}}\times 0.00346\times 0.698=0.03\text{\hspace{0.17em}}\text{m}/\text{sec}$

On the other hand, for the second and easier route, imagine increasing the space symmetrically. The flow is transformed into laminar flow in a duct with a
$gap=\stackrel{\xaf}{\alpha \xi}=2\times 0.0208\text{\hspace{0.17em}}\text{m}$ and a velocity distribution *α*, *β*, *ξ*.

The flow now is symmetric and so,

$\sigma =\left|\begin{array}{cc}0& 0.00208\\ 0.00208& 0\end{array}\right|$

And back to Equation (27a),

${V}^{\u2033}=\frac{2\times \left(2\times 0.0208\right)}{8\times 0.5\times {10}^{-3}}\times 0.00208\times 0.69=0.03\text{\hspace{0.17em}}\text{m}/\text{sec}$

Finally, both methods provide the same result.

Crosschecking this value with that from the exact analytic solution [7],

${V}^{\u2033}=\frac{\nabla p}{12{\mu}^{\u2033}}{b}^{2}\left(\frac{{\mu}^{\prime}+7{\mu}^{\u2033}}{{\mu}^{\prime}+{\mu}^{\u2033}}\right)=0.031\text{\hspace{0.17em}}\text{m}/\text{sec}$

Find the average velocity ${V}^{\prime}$ in the section where the fluid has viscosity ${\mu}^{\prime}$.

From Figure 7, appears that in the cases where there is a substantial difference in the fluid viscosities, *i*.*e*.,
${\mu}^{\prime}\gg {\mu}^{\u2033}$, the velocity gradient at the point *γ* is almost zero. This means that one can imagine a symmetric laminar flow in a duct with the total gap equal to 2*b*. The tensor in this case is symmetric,

$\sigma =\left|\begin{array}{cc}0& {{\tau}^{\u2033}}_{w}\\ {{\tau}^{\u2033}}_{w}& 0\end{array}\right|$ *and
${{\tau}^{\u2033}}_{w}=\nabla p\left(b+L\right)=0.00292$ *

So, again from Equation (27a),

${V}^{\prime}=\frac{2\times \left(2\times 0.025\right)}{8\times {10}^{-3}}\times 0.00292\times 0.698=0.0255\text{\hspace{0.17em}}\text{m}/\text{sec}$

And the exact analytic solution gives [7],

${V}^{\prime}=\frac{\nabla p}{12{\mu}^{\prime}}{b}^{2}\left(\frac{7{\mu}^{\prime}+{\mu}^{\u2033}}{{\mu}^{\prime}+{\mu}^{\u2033}}\right)=0.026\text{\hspace{0.17em}}\text{m}/\text{sec}$

Moreover, given that the calculations are based on symmetric laminar flow in a duct and, hence, the velocity distribution is almost parabolic, it is possible to calculate even the maximum velocity as simply as,

${{V}^{\u2033}}_{\mathrm{max}}=\frac{3}{2}{V}^{\u2033}=\frac{3}{2}\times 0.03=0.045\text{\hspace{0.17em}}\text{m}/\text{sec}$

And from the analytical solution [7],

${{V}^{\u2033}}_{\mathrm{max}}=\frac{\nabla p}{2{\mu}^{\u2033}}{b}^{2}\left(\frac{2{\mu}^{\u2033}}{{\mu}^{\prime}+{\mu}^{\u2033}}+\frac{{L}^{2}}{{b}^{2}}\right)=0.0434\text{\hspace{0.17em}}\text{m}/\text{sec}$

Respectively, for*
${{V}^{\prime}}_{\mathrm{max}}$ *.

Moreover, given velocity’s distribution is almost a parabola of which the apex ( ${{V}^{\u2033}}_{\mathrm{max}}$ ) and the base are known it is feasible to build back the parabola and calculate the velocity at any point with good approximation.

*Unfortunately*,*a check of this method*,*for a turbulent flow and at least for the part
${\mu}^{\u2033}$ * *has not been made due to the lack of experimental data*.

5. Model Validation

In doing so the proposed model (practically Equation (45)) has been checked against the method proposed by Jones *et al.* [4] and generally accepted [2] and with some experimental data [4] [10]. Obviously for the laminar flow the check was realized against the analytically obtained exact solutions.

Suppose there is a situation to calculate pressure gradient in an annulus. The geometry of the annulus is:
$K=0.8$ annular flow of water*
$R{e}_{h}={10}^{5}$ *,
$\rho =1000\text{\hspace{0.17em}}\text{kg}/{\text{m}}^{\text{3}}$,
$\mu ={10}^{-3}\text{kg}/\text{m}\cdot \mathrm{sec}$,*
$\lambda =0.9$ *,*
$R=0.05\text{\hspace{0.17em}}\text{m}$ *, *hence*
$V=5\text{\hspace{0.17em}}\text{m}/\text{sec}$.

1) Jones (*et al.*) approach. For
$K=0.8$ they propose a correction factor*
${\Phi}^{*}=0.667$ *, (*see *Appendix B Equation (B.7)).Hence,
$R{e}_{{d}_{h}}\times {\Phi}^{*}=R{e}_{{d}_{h}}^{\ast}$,
$R{e}_{{d}_{h}}^{\ast}=66700$. With this number one can enter to the Moody’s chart to find the friction factor. Alternatively, there is a chart prepared with experimental data by Leung [10] for an annulus with*
$K=0.8$ * (also in Figure 8 in [4] ) and can be used.

Both give, $f=0.019$ and from Equation (39),

$5=7.25\times \sqrt{\frac{2\times 0.02}{1000}}\sqrt{\nabla p}$ hence $\nabla p=11980\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$.

2) Implementing the present method and $R{e}_{h}={10}^{5}$ enter Moody’s or use Equation (40). Find friction factor $f=0.0178$. Enter directly to Equation (45) and find the gradient.

$5=2\times 7.49\times 0.00317\times 0.974\times \sqrt{\nabla p}$ hence $\nabla p=11686\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$

Hence, this result identifies with the previous based on real data [10] and/or the well accepted Jones (*et al.*) method.

The present approach is more straightforward, no need to define and calculate the correction factor
${\Phi}^{*}$ or prepare a chart for every different *Κ*.

Repeating as above but for $R{e}_{h}={10}^{4}$ the results are respectively,

$\nabla p=208\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$ and $\nabla p=205\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$.

6. Conclusions and Discussion

It has been shown, herein, that in a non-symmetrical flow just the use of the Hydraulic Diameter concept is not sufficient to describe the flow. More information is required. The approach is based on classical mechanics, while the two flow types, laminar and turbulent, have been considered separately.

In laminar flow, the fluid was considered as a continuum with cohesive forces present. So, the stresses, instead of deforming an element of the fluid, can move it. With the help of Eigenvectors-Eigenvalues, a direction was found where shear stresses disappear. Finally, the driving force has been calculated. It turned out that the net moving force forward is just the Z component of the driving force, *i*.*e*., to fully describe the flow, the value of the angle
$\theta $ is needed. Finally, Equation (27) and/or Equation (28) relate the average velocity to the pressure gradient in a straightforward way.

In turbulent flow, as it has been explained in detail, there are no cohesive forces and therefore, the previous mechanical approach has no sense. Fortunately, the Eigenvector keeps indicating the point of maximum velocity. This point is a function of
$\lambda $. *I*.*e*., the point of maximum velocity or zero shear stresses is just a function of the geometry of the annuli. This very fact also explains the success of the “Laminar equivalent diameter” concept. Both laminar and turbulent, although completely different flows, share the same annuli geometry.

Another possible explanation could be that the present approach is focused on the wall stresses and more precisely to their difference. Nevertheless, a different but real stress distribution still exists in turbulent flow. Hence, a point of maximum velocity does exist in the turbulent case as well.

To point out that in this approach the stress tensor has been calculated on the non-symmetrical wall stresses. For more complicated cases, the stress tensor can be defined at any point in the core of the liquid. This tensor of course will be function of the point, exactly like the metric tensor on curved surface in the differential geometry.

Definitely Equation (45) relates with great accuracy the average velocity to the pressure gradient in turbulent flow and for any kind of concentric annuli. At the same time, the eigenvectors show the direction of the driving force in the core of the fluid.

The advantage against the existing method to calculate average velocity in annuli is that there is no need to deal with correction coefficients etc. or “equivalent hydraulic diameter” anymore. In fact, it is positive to get rid of definitions that complicate with no reason the already difficult fluid mechanics science.

Hence, it is had been shown that just the use of the “Hydraulic Diameter” concept is not sufficient for non-symmetric flows like in annuli.

Instead, was found that for either laminar or turbulent flow in rectangular or square ducts, the concept of “Hydraulic Diameter” is sufficient. This fact is in line with the topic of the present paper given the fact that the flow in a duct is symmetric. Hence, several new equations have been provided. Equation (46) is valid for any rectangular duct in turbulent flow. Equation (47) is valid for any square duct in turbulent flow as well. Moreover, both Equation (27a) and Equation (27b) are valid for any rectangular duct in laminar flow.

Therefore, again has been shown that there is no need of the suggested route [2] [3] along with the use of “laminar equivalent diameter” concept and further complications for the rectangular ducts.

To conclude the approach, based on the “laminar equivalent diameter” concept proposed by Jones *et al.* [3] [4] and generally accepted [2], is more complicated, for no reason than the present one and with a non-robust theoretical argument.

Moreover, the present original method, based on Eigenvalues and Eigenvectors, is a promising method to simplify problems in fluid mechanics. In addition, it will help engineers to comprehend fluid dynamics in a different way. This could be very useful and, although most of these problems are already successfully solved, a deeper understanding is always welcome. Additionally, heat transfer, mass transfer and others are strongly linked to fluid dynamics.

In addition, the experts in the field could possibly benefit from the present approach in dealing with non-symmetric flow like, flow of immiscible liquids, flow in a curved conduit, flow around non-symmetric objects, flow with non-symmetric heat exchange, flow with chemical reactions etc. [11] [12] [13] [14].

Nomenclature

*A* = cross section area (m^{2})

*d _{h}* or

*f* = friction factor

*h* = half a duct’s gap (m)

*K* = radius ratio *R _{i}*/

*L* = length (m)

*ʟ* = Max. velocity coordinate to adjacent immiscible fluids flow (m)

*P* = pressure (N/m^{2})

*P* = cross-section wetted perimeter (m)

*r* = the radial coordinate (m)

*R* = the radius (m)

*Re* = Reynolds’ number

*Re _{h}* =Reynolds’ number defined using

*Re ^{*}* = modified Reynolds’ number

*V* = average velocity (m/sec)

${v}_{1},{v}_{2}$ = the eigenvectors

*ε*/*d* = roughness ratio

${\lambda}_{1},{\lambda}_{2}$ = the eigenvalues

$\lambda $ (in bold) defines as $\lambda R$ the coordinate of max. velocity

*μ* = dynamic viscosity (kg/msec)

*ρ* = density (kg/m^{3})

*σ* = stress (N/m^{2})

*τ* = shear stress (N/m^{2})

${\Phi}^{*}$ = geometry coefficient in [2] and [3]

Appendix A

Calculating $\lambda $

Imagine there is need to transform a biplane into a single wing plane (see FigureA1).

In doing so, it is necessary not just the area of the new wing to be equivalent to the two ex-areas, but also there is need to consider the position and the distance from the center of mass.

Same reasoning applies in the concentric tubes and thus,

$\left(2\pi \lambda RL\right)\lambda R=\frac{\left(2\pi RL\right)R-\left(2\pi KRL\right)KR}{\mathrm{ln}\frac{\left(2\pi RL\right)R}{\left(2\pi KRL\right)KR}}$ $\to $ ${\lambda}^{2}=\frac{1-{K}^{2}}{\mathrm{ln}\frac{1}{{K}^{2}}}$

*I*.*e*., *the* *wind* *passes* *through* *the* *wings* *and* *creates* *moments* *relative* *to* *the* *mass* *center* *and* *proportional* *to* *the* *areas* *involved*. *There* *is* *a* *similar* *situation* *in* *the* *annulus*. *Practically* *there* *is* *need* *to* *calculate* *an* *equivalent* *or* *an* *average* *of* *the* *first* *area* *moment*. *In* *fact*, *the* “*logarithmic* *mean*” *is* *the* *most* *appropriate* *to* *be* *used* *in* *such* *geometry* *problems*.

Figure A1. Transforming a biplane into a single wing plane.

Appendix B

Trying to explain the adopted method in [4].

In paper [4], Darcy’s definition for the friction factor is used,

$\frac{1}{2}\rho {V}^{2}=\frac{4{\tau}_{w}}{f}=\frac{4}{f}\frac{\nabla p}{4}D$

${V}^{2}=\frac{1}{f}\frac{2D}{\rho}\nabla p$ (B.1)

Together with the semiempirical equation for the laminar flow,

$f=\frac{64}{Re}$ (B.2)

Accepting the definition of ${D}_{h}$ in annuli,

${D}_{h}=2\left(R-KR\right)=2R\left(1-K\right)$ (B.3)

From Equation (B.1), Equation (B.2) and Equation (B.3) we have,

$V=\frac{2{D}^{2}{\left(1-K\right)}^{2}}{64\mu}\nabla p$ (B.4)

Considering the analytical solution as well [7],

$V=\frac{\nabla p}{8\mu}{R}^{2}\left(\frac{1-{K}^{4}}{1-{K}^{2}}-\frac{1-{K}^{2}}{\mathrm{ln}\frac{1}{K}}\right)$ (B.5)

Or

$V=\frac{\nabla p}{8\mu}{R}^{2}\left(1+{K}^{2}-\frac{1-{K}^{2}}{\mathrm{ln}\frac{1}{K}}\right)$ (B.6)

Everything in the parenthesis has to do with geometry.

In [4], a correction factor is defined as,

${\Phi}^{*}=\frac{1}{{\left(1-K\right)}^{2}}\left(1+{K}^{2}-\frac{1-{K}^{2}}{\mathrm{ln}\frac{1}{K}}\right)$ (B.7)

Now multiplying the semiempirical Equation (B4) with Equation (B7) we get,

$V=\frac{2{D}^{2}\left(1-{K}^{2}\right)}{64\mu}\nabla p{\Phi}^{*}=\frac{\nabla p}{8\mu}{R}^{2}\left(1+{K}^{2}-\frac{1-{K}^{2}}{\mathrm{ln}\frac{1}{K}}\right)$

*I*.*e*., *with* *this* *definition* *of* *the* *correction* *factor*
${\Phi}^{*}$ *it* *has* *been* *possible* *to* *transform* *a* *semiempirical* *equation* *to* *an* *analytical* *one* *at* *least* *for* *the* *laminar* *flow*. *Extrapolating* *this* *idea* *to* *the* *turbulent* *flow* *seems* *working* *well* *enough*. *Keeping* *in* *mind* *that* Equation (B.2) *is* *based* *on* *experiments* *to* *a* *circular* *cross* *section* *pipe*.

The corresponding correction factor by Jones for flow in a rectangular duct [3] is given by,

${\Phi}^{*}~\frac{2}{3}+\frac{11}{24}\frac{s}{w}\left(2-\frac{s}{w}\right)$ (B.8)

Appendix C

Herein, the derivation and the validity of Equation (28) will be explained.

Consider the case of $K=0.2$ water in laminar flow

$R=0.05\text{\hspace{0.17em}}\text{m}$, $\lambda =0.546$, $\nabla p=1\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$

${\tau}_{o}^{w}=\frac{\nabla p}{2}R\left(1-{\lambda}^{2}\right)=\nabla p\times 0.0175\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

${\tau}_{i}^{w}=\frac{\nabla p}{2}R\left(K-\frac{{\lambda}^{2}}{K}\right)=\nabla p\times 0.0323\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

$\sigma =\left|\begin{array}{cc}0& {\tau}_{i}^{w}\\ {\tau}_{o}^{w}& 0\end{array}\right|=\nabla p\left|\begin{array}{cc}0& 0.0323\\ 0.0175& 0\end{array}\right|$

${\sigma}_{\text{diagonal}}=\nabla p\left|\begin{array}{cc}{\lambda}_{1}& 0\\ 0& {\lambda}_{2}\end{array}\right|=\nabla p\left|\begin{array}{cc}0.0238& 0\\ 0& -0,0238\end{array}\right|$

${v}_{1}=\left|\begin{array}{c}1.3586\\ 1\end{array}\right|$ and ${v}_{2}=\left|\begin{array}{c}-1.3586\\ 1\end{array}\right|$

And $\theta ={53.6}^{\circ}$

Using Equation (28) to calculate the average velocity,

$V=\frac{{D}_{h}}{8\mu}\lambda \mathrm{cos}\theta \text{\hspace{0.05em}}\nabla p$ (28)

$V=\frac{0.08\times 0.0238\times 0.593}{8\times {10}^{-3}}\nabla p$

$V=0.141\text{\hspace{0.17em}}\text{m}/\text{sec}$

Excellent fit with the analytical solution.

It is necessary to try to understand what is happening.

From FigureC1, we see that the angle
$\theta ={53.6}^{\circ}$ indicates the direction to point *M*, which is the coordinate of maximum velocity and/or zero stresses.

$T=\left|\begin{array}{c}{T}_{r}\\ {T}_{z}\end{array}\right|=0.0238\left|\begin{array}{c}\mathrm{sin}\theta \\ \mathrm{cos}\theta \end{array}\right|=0.0238\left|\begin{array}{c}0.805\\ 0.593\end{array}\right|$

$\left|\begin{array}{c}{T}_{r}\\ {T}_{z}\end{array}\right|=\left|\begin{array}{c}0.019\\ 0.0141\end{array}\right|$ $\to $ ${T}_{z}=0.0141\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

It seems that this traction vector component is an average value of what is happening in the fluid.

Finding the angle *φ* of the two normalized eigenvectors,

$\mathrm{cos}\phi =\frac{{v}_{1}\circ {v}_{2}}{1\times 1}$

${v}_{1}\circ {v}_{2}=-{0.805}^{2}+{0.593}^{2}=-0.2964$

$\phi ={107.24}^{\circ}$

Therefore, the angle of the two eigenvectors is 107.24˚. So, the eigenvectors are not orthogonal, which is reasonable given that the stress tensor is not symmetric. Given the above, the construction of a new triangle (*AM*’*E*) *is* *needed* (*see* FigureC2). For this new triangle the angle, the base, and the point *N* are known.

The triangle (*AM*’*E*) can be constructed graphically using a protractor or by solving the following implicit equation,

$\frac{\mathrm{tan}\measuredangle A}{\mathrm{tan}\measuredangle E}=\frac{0.0225}{0.0175}=1.286$

along with,

$\measuredangle A+\measuredangle E={72.76}^{\circ}$

Figure C1. Annular flow, K = 0.2, R = 0.05 m. Eigenvectors supposed to be orthogonal.

Figure C2. Annular flow, K = 0.2, R = 0.05 m. Eigenvectors not orthogonal.

Figure C3. Annular flow, K = 0.2, R = 0.05 m. Traction vectors not orthogonal.

Solving the system by trial and error,

$\measuredangle A={39.5}^{\circ}$, $\measuredangle E={32.5}^{\circ}$, $N{M}^{\prime}=0.0\text{143}\text{\hspace{0.17em}}\text{m}$

In FigureC3, the two eigenvectors are aligned with the traction vectors ${T}_{1},{T}_{2}$. Given that the eigenvectors are normal (in fact they define it) to the principal planes, the traction vectors are also normal to the principal planes.

${T}_{1}=\left|\begin{array}{c}{T}_{1r}\\ {T}_{1z}\end{array}\right|=\sigma \cdot {v}_{1}=\lambda \left|\begin{array}{c}\mathrm{sin}{57.5}^{\circ}\\ \mathrm{cos}{57.5}^{\circ}\end{array}\right|$

or

${T}_{1}=\left|\begin{array}{c}{T}_{1r}\\ {T}_{1z}\end{array}\right|=0.0238\left|\begin{array}{c}0.843\\ 0.537\end{array}\right|$

${T}_{1z}=0.0238\times 0.537=0.0128\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

In the same way,

${T}_{2z}=0.0238\times 0.636=0.0151\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

And the geometric mean,

${T}_{aver.z}=\sqrt{0.0128\times 0.0151}=0.014\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

So, the average value of the Z component of the two traction vectors identifies with the value obtained from Equation (28).

In fact, the angle $\theta ={53.6}^{\circ}$ is half of $\phi ={107.24}^{\circ}$

Moreover, the average of the complementary angles to
$\measuredangle A,\measuredangle E$ is exactly 53.6˚. *I*.*e*., the two non-orthogonal eigenvectors define the point *M*’ into the 2D space.

Appendix D

What happens in the cas*e* of
$K\to 0$ ?

From Equation (22), $\lambda \to 0$ and $\theta ={90}^{\circ}$

From FigureD1, now ${\tau}_{i}=0$

And the stress tensor must be,

$\sigma =\frac{\nabla p}{2}R\left|\text{\hspace{0.05em}}\begin{array}{cc}0& 0\\ 1-{\lambda}^{2}& 0\end{array}\text{\hspace{0.05em}}\right|$

And for $\lambda =0$

$\sigma =\frac{\nabla p}{2}R\left|\text{\hspace{0.05em}}\begin{array}{cc}0& 0\\ 1& 0\end{array}\text{\hspace{0.05em}}\right|$

This tensor is not diagonalizable, but without changing the substance of the tensor, filling back with the static pressure assuming to be 1 N/m^{2}^{ }

Figure D1. Flow in a circular tube with $R=0.05\text{\hspace{0.17em}}\text{m}$.

$\sigma =\frac{\nabla p}{2}R\left|\begin{array}{cc}-1& 0\\ 1& -1\end{array}\right|$

This tensor is still not diagonalizable, but has a non-zero determinant, one eigenvalue $\lambda =-1$ and just one eigenvector, ${v}_{1}=\left|\text{\hspace{0.05em}}\begin{array}{c}0\\ 1\end{array}\text{\hspace{0.05em}}\right|$.

So, $\theta =0$ and $\mathrm{cos}\theta =1$.

From Equation (27),

$V=\frac{{D}_{h}}{8\mu}\sqrt{\left|\Vert \sigma \Vert \right|}\mathrm{cos}\theta $ $\to $ $V=\frac{{R}^{2}}{8\mu}\nabla p$

*I*.*e*., *the* *analytic* *equation* *for* *laminar* *flow* *in* *a* *circular* *tube*.

What about turbulent flow? Going back to the definition in Equation (31),

$\frac{1}{2}\rho {V}^{2}=\frac{1}{f}\sqrt{\left|\Vert \sigma \Vert \right|}$

$\frac{1}{2}\rho {V}^{2}=\frac{1}{f}\frac{\nabla p}{2}R$ $\to $ $V=\frac{1}{\sqrt{f}}\sqrt{\frac{R}{\rho}}\sqrt{\nabla p}$

Inserting coefficient 2 to fix the problem with the friction coefficient definition,

$V=2\frac{1}{\sqrt{f}}\sqrt{\frac{R}{\rho}}\sqrt{\nabla p}$ or $V=\frac{1}{\sqrt{f}}\sqrt{\frac{2D}{\rho}}\sqrt{\nabla p}$

This last equation identifies with Equation (39); the classic one for turbulent flow in a circular pipe.

Appendix E

*Solving* *the* *problem* (2*A-*3) *page* 62 [7].

To better understand the use of Equation (27) for laminar flow in an annulus, the above problem is solved. The units are converted to SI.

Data: ${R}_{i}=0.0124\text{\hspace{0.17em}}\text{m}$, ${R}_{o}=0.0275\text{\hspace{0.17em}}\text{m}$

$\rho =80.3\text{\hspace{0.17em}}\text{lb}/{\text{ft}}^{\text{3}}$, $\mu =0.05738\left(\text{kg}/\text{m}\cdot \mathrm{sec}\right)$

$\Delta p=37163\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$, $L=8.23\text{\hspace{0.17em}}\text{m}$

Find: $Q=?\text{\hspace{0.17em}}{\text{m}}^{\text{3}}/\text{sec}$

So, $K=\frac{{R}_{i}}{{R}_{o}}=0.45$

From Equation (2), $\lambda =0.706$

$\nabla p=\frac{\Delta p}{L}=4516\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{3}}$

From Equation (5) and Equation (6),

${\tau}_{i}^{w}=41\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

${\tau}_{o}^{w}=31.1\text{\hspace{0.17em}}\text{N}/{\text{m}}^{\text{2}}$

Ready to build the tensor,

$\sigma =\left|\begin{array}{cc}0& {\tau}_{i}^{w}\\ {\tau}_{o}^{w}& 0\end{array}\right|=\left|\begin{array}{cc}0& 41\\ 31.1& 0\end{array}\right|$

There are many online algorithms to diagonalize the tensor $\sigma $,

${\sigma}_{\text{diagonal}}=\left|\begin{array}{cc}{\lambda}_{1}& 0\\ 0& {\lambda}_{2}\end{array}\right|=\left|\begin{array}{cc}35.7& 0\\ 0& -35.7\end{array}\right|$

And the eigenvectors are,

${v}_{1}=\left|\begin{array}{c}1.148\\ 1\end{array}\right|$ and ${v}_{2}=\left|\begin{array}{c}-1.148\\ 1\end{array}\right|$

*
${\mathrm{tan}}^{-1}\theta =1.148$ *
$\to $
$\theta ={48.9}^{\circ}$

Or, $\mathrm{cos}\theta =0.657$

$V=\frac{{D}_{h}}{8\mu}\sqrt{\left|\Vert \sigma \Vert \right|}\mathrm{cos}\theta $ (27)

For an annulus, the Hydraulic Diameter is 2 times the gap,

${D}_{h}=2R\left(1-K\right)=0.03\text{\hspace{0.17em}}\text{m}$

The determinant of
${\sigma}_{\text{diagonal}}$ or
$\sigma $ *is*,

$\Vert \sigma \Vert =35.7\times \left(-35.7\right)$

So,

*
$V=\frac{0.03}{8\times 0.05738}\times 35.7\times 0.657=1.53\text{\hspace{0.17em}}\text{m}/\text{sec}$ *

*
$Q=1.53\times A=1.53\times 0.00189=0.003\left({\text{m}}^{\text{3}}/\text{sec}\right)$ *

Or $Q=0.11\left({\text{ft}}^{\text{3}}/\mathrm{sec}\right)$.

In the case a computer algorithm is not available, how one can find the ${\sigma}_{\text{diagonal}}$ and $\theta $ ?

Observing that,

$35.7=\sqrt{41\times 31.1}$

*I*.*e*., the quantity 35.7 is just the geometric average of the two shear stresses.

To find the angle $\theta $, the classical formula which provides the rotation angle to the principal planes in the “Strength of Materials” theory, is used.

$\mathrm{tan}2\theta =\frac{{\tau}_{12}}{\frac{{\sigma}_{1}-{\sigma}_{2}}{2}}$ (E.1)

Equation (E.1) is useful for symmetric tensors for which, ${\tau}_{12}={\tau}_{21}$.

In this case the tensor is not symmetric, however, one can imagine he is already on the principal plains, and wants to rotate back,

$\mathrm{tan}2\left(90-\theta \right)=\frac{\sigma}{\frac{{\tau}_{12}-{\tau}_{21}}{2}}$ (E.2)

Hence, Equation (E.2) becomes,

$\mathrm{tan}2\left(90-\theta \right)=\frac{35.7}{\frac{41-31,1}{2}}$ $\to $ $\theta ={48.9}^{\circ}$

Alternatively, in cases $K\ge 0.3$, use Equation (22).

Finally, now a formula based on continuum mechanics can be introduced that calculates the average velocity in any annulus in laminar flow,

$V=\frac{{R}^{2}\left(1-K\right)}{8\mu}\nabla p\sqrt{\left(1-{\lambda}^{2}\right)\left(K-\frac{{\lambda}^{2}}{K}\right)}\mathrm{cos}\theta $ (E.3)

Along with,

$\mathrm{tan}2\left(90-\theta \right)=2\frac{\sqrt{\left(1-{\lambda}^{2}\right)\left(K-\frac{{\lambda}^{2}}{K}\right)}}{\left(K-\frac{{\lambda}^{2}}{K}\right)-\left(1-{\lambda}^{2}\right)}$ (E.4)

Conflicts of Interest

The authors declare no conflicts of interest regarding the publication of this paper.

[1] |
Colebrook, C.R. (1939) Turbulent Flow in Pipes with Particular Reference to the Transition Region between the Smooth and the Rough Pipe Laws. Journal of the Institution of Civil Engineers (London), 11, 133-156. https://doi.org/10.1680/ijoti.1939.13150 |

[2] | White, M.F. (2015) Fluid Mechanics. 8th Edition in SI Units, First Published Dec. |

[3] |
Jones, Jr., O.C. (1976) An Improvement in the Calculation of Turbulent Friction in Rectangular Ducts. ASME Journal of Fluids Engineering, 98, 173-181. https://doi.org/10.1115/1.3448250 |

[4] |
Jones, Jr., O.C. and Leung, J.C.M. (1981) An Improvement in the Calculation of Turbulent Friction in Smooth Concentric Annuli. ASME Journal of Fluids Engineering, 103, 615-623. https://doi.org/10.1115/1.3241781 |

[5] | Kratz, A.P., Macintire, H.J. and Gould, R.E. (1930) Flow of Liquids in Pipes of Circular and Annular Cross-Sections. Bulletin No. 222, University of Illinois, Champaign, IL. |

[6] |
Caetano, E.F., Shoham, O. and Brill, J.P. (1992) Upward Vertical Two-Phase Flow through an Annulus-Part I: Single-phase Friction Factor, Taylor Bubble Rise Velocity and Flow Pattern Prediction. The Journal of Energy Resources Technology, 114, 1-13. https://doi.org/10.1115/1.2905917 |

[7] | Bird, R.B., Stewart, W.E. and Lightfoot, E.L. (2006) Transport Phenomena. 2nd Edition, John Wiley & Sons, Inc., New York. |

[8] |
Haaland, S.E. (1983) Simple and Explicit Formulas for the Friction Factor in Turbulent Pipe Flow. The Journal of Fluids Engineering, 105, 89-90. https://doi.org/10.1115/1.3240948 |

[9] |
Wang, C.Y. (2017) Starting Flow in a Channel with two Immiscible Fluids. The Journal of Fluids Engineering, 139, Article ID: 124501. https://doi.org/10.1115/1.4037495 |

[10] |
Leung, J.C.M. (1975) Occurrence of Critical Heat Flux during Blowdown with Flow Reversal. M.S. Thesis, Northwestern University, Evanston, IL. https://doi.org/10.2172/7246823 |

[11] |
Tanveer, A., Hayat, T. and Alsaedi, A. (2021) Numerical Simulation for Peristalsis of Sisko Nanofluid in Curved Channel with Double-Diffusive Convection. Ain Shams Engineering Journal, 12, No. 3. https://doi.org/10.1016/j.asej.2020.12.019 |

[12] |
Tanveer, A., Salahuddin, T., Khan, M., Malik, M.Y. and Alqarni, M.S. (2020) Theoretical Analysis of Non-Newtonian Blood Flow in a Microchannel. Computer Methods and Programs in Biomedicine, 191, Article ID: 105280. https://doi.org/10.1016/j.cmpb.2019.105280 |

[13] |
Tanveer, A., Khan, M., Salahuddin, T., Malik, M.Y. and Khan, F. (2019) Theoretical Investigation of Peristaltic Activity in MHD Based Blood Flow of Non-Newtonian Material. Computer Methods and Programs in Biomedicine, 187, Article ID: 105225. https://doi.org/10.1016/j.cmpb.2019.105225 |

[14] |
Khan, M., Shahid, A., El Shafey, M., Salahuddin, T. and Khan, F. (2020) Predicting Entropy Generation in Flow of Non-Newtonian Flow due to a Stretching Sheet with Chemically Reactive Species. Computer Methods and Programs in Biomedicine, 181, Article ID: 105246. https://doi.org/10.1016/j.cmpb.2019.105246 |

Journals Menu

Contact us

+1 323-425-8868 | |

customer@scirp.org | |

+86 18163351462(WhatsApp) | |

1655362766 | |

Paper Publishing WeChat |

Copyright © 2024 by authors and Scientific Research Publishing Inc.

This work and the related PDF file are licensed under a Creative Commons Attribution 4.0 International License.