The problem
A Bingham plastic fluid flows upwards axially through an inclined circular pipe of length $L$ and radii $R$. The pressures at $z=0,L$ are $p=P_0,P_L$, respectively. Calculate the velocity profile in steady state and the corresponding volumetric flow rate. Exit and entrances losses, and heating due to friction are very small, so that these can be omitted.
 |
| Fig. 1 Problem sketch for the Bingham plastic plug flow through an inclined pipe |
The solution for the velocity distribution and flow rate of a Bingham plastic fluid
According to the sketch in Fig. 1, the geometry of the problem is cylindrical. So, it seems obvious to use cylindrical coordinates for the solution of the problem.
Assumptions and physics of the fluid flow of a Bingham plastic
Next, according to the physics of the problem the following conclusions are made,
- the fluid moves only along the axial direction from $z=0$ to $z=L$. This is, the fluid moves parallel to the cylinder axis;
- the fluid moves upwards through the inclined pipe;
- the fluid does not move radial or angularly. From a transversal view, the fluid does not move left, right or in another radial direction. Also, the fluid does not circle motions;
- considering that the function representing the fluid velocity is $\textbf{v}=(v_r,v_\theta, v_z)$ and in agreement with the above statements $v_r=0$, $v_\theta=0$ and $v_z\ne 0$ follow;
- the component $v_z$ does not vanish and may depend on $(r,\theta,z,t)$. In this case, it only depends on $r$. One arrives to this conclusion as follows. The continuity equation assures that the fluid velocity at $z=0$ and $z=L$ is the same, so that it does not depends on $z$. Besides, it is easy to realize that the fluid velocity remains unchanged under angular variations;
- stationary conditions are considered. This means that there is no acceleration and that the derivatives of $\textbf{v}$ are zero or that there is no temporal dependency;
- the pressure difference between $z=0$ and $z=L$ sets the fluid in motion;
- the pressure $p$ is in general a function of $(r, \theta, z,t)$; but in this case it only depends on $z$. One arrives to this conclusion as follows. Since the fluid is not moving radial nor angularly, then there is no dependency of the pressure on $(r,\theta)$;
- the pressure gradient is generated externally and not shown in Fig. 1. However, since the fluid moves upwards the acceleration due to gravity has the opposite effect;
- the problem has angular symmetry. This is, none of the variables of the problem as the pressure, the velocity or the stress depend on $\theta$. In other words, it does not matter the assigned value to $\theta$, the system will not change. $\theta=0$ will be used;
- the effect of temperature on the fluid properties is very small and will be discarded. Then, $\rho$ and $\mu$ are constants;
- since this is a Bingham plastic fluid it will only move if a critical stress, sometimes named \textit{yield stress}, is surpassed;
- due to the viscosity and to the plastic nature of the fluid it is expected that the fluid in the pipe experiences great deformation near the wall; but little deformation around the center area of the pipe (see Fig. (2) for more details);
- if the stress applied to the fluid, to set it in motion, is too high the fluid may behave as a newtonian fluid. In other words, it would be the case of $\mathbf{\tau}>>\mathbf{\tau_0}$.
 |
| Fig. 2 A representation of how the plastic Bingham fluid is deformed when moving through a pipe |
In order to determine the velocity distribution two equations should be considered: one defining the type of fluid and another one that models the fluid motion under any circumstance. These equations are,
- the fluid constitutive equation.- Since the working fluid is a plastic of Bingham the corresponding equation is,
$\text{if}\quad \mathbf{\tau}\geq\mathbf{\tau_0}, \quad \mathbf{\tau}=\mathbf{\tau_0} - \mu_\infty \left[\nabla \mathbf{v} + \left(\nabla \mathbf{v}\right)^T\right]$ Eq. (01)
The constitutive Eq. (01) is different from that for newtonian and power law fluids due to its split form. Although it looks complex it is a very simple way to represent the clay-like behavior of a material commonly found in industry. A rough pictorial representation of Eq. (01) is shown in Fig. 2 where the plug flow corresponds to the zero strain or deformation part of Eq. (01), when $\mathbf{\tau}\leq\mathbf{\tau_0}$, while the parabolic flow section corresponds to the stress-strain relationship given in Eq. (01) when $\mathbf{\tau}\geq\mathbf{\tau_0}$.
Next, Eq. (01) is rewritten in extended form using cylindrical coordinates as,
2\dfrac{\partial v_r}{\partial r} & r\dfrac{\partial}{\partial r}\left(\dfrac{v_\theta}{r}\right)+\dfrac{1}{r}\dfrac{\partial v_r}{\partial \theta} & \dfrac{\partial v_r}{\partial z}+\dfrac{\partial v_z}{\partial r} \\
r\dfrac{\partial}{\partial r}\left(\dfrac{v_\theta}{r}\right)+\dfrac{1}{r}\dfrac{\partial v_r}{\partial \theta} & 2\left(\dfrac{1}{r}\dfrac{\partial v_\theta}{\partial \theta}+\dfrac{v_r}{r}\right) & \dfrac{1}{r}\dfrac{\partial v_z}{\partial \theta}+\dfrac{\partial v_\theta}{\partial z} \\
\dfrac{\partial v_r}{\partial z}+\dfrac{\partial v_z}{\partial r} & \dfrac{1}{r}\dfrac{\partial v_z}{\partial \theta}+\dfrac{\partial v_\theta}{\partial z} & 2\dfrac{\partial v_z}{\partial z}
\end{bmatrix} =0$ Eq. (02)
for the case of $\mathbf{\tau}\leq\mathbf{\tau_0}$ where no deformation occurs. For the case of $\mathbf{\tau}\geq\mathbf{\tau_0}$ the equation becomes,
\tau_{rr}\cos^2 \theta-2\tau_{r\theta}\sin \theta \cos \theta+\tau_{\theta \theta}\sin^2 \theta & \tau_{rr}\sin \theta \cos \theta + \tau_{r\theta}\left[\cos^2 \theta-\sin^2 \theta\right]-\tau_{\theta \theta}\sin \theta\cos \theta & \tau_{rz}\cos \theta -\tau_{z\theta}\sin \theta\\
\tau_{rr}\sin \theta \cos \theta+\tau_{r\theta}\left[\cos^2 \theta-\sin^2 \theta\right]-\tau_{\theta \theta}\sin \theta\cos \theta & \tau_{rr}\sin^2 \theta+2\tau_{r\theta}\sin \theta \cos \theta+\tau_{\theta \theta}\cos^2 \theta & \tau_{rz}\sin \theta +\tau_{z\theta}\cos \theta\\
\tau_{rz}\cos \theta -\tau_{z\theta}\sin \theta & \tau_{rz}\sin \theta +\tau_{z\theta}\cos \theta & \tau_{zz}
\end{bmatrix}\nonumber\\
=\begin{bmatrix}
\tau_0 & \tau_0 & \tau_0 \\
\tau_0 & \tau_0 & \tau_0 \\
\tau_0 & \tau_0 & \tau_0
\end{bmatrix}
-\mu_\infty
\begin{bmatrix}
2\dfrac{\partial v_r}{\partial r} & r\dfrac{\partial}{\partial r}\left(\dfrac{v_\theta}{r}\right)+\dfrac{1}{r}\dfrac{\partial v_r}{\partial \theta} & \dfrac{\partial v_r}{\partial z}+\dfrac{\partial v_z}{\partial r} \\
r\dfrac{\partial}{\partial r}\left(\dfrac{v_\theta}{r}\right)+\dfrac{1}{r}\dfrac{\partial v_r}{\partial \theta} & 2\left(\dfrac{1}{r}\dfrac{\partial v_\theta}{\partial \theta}+\dfrac{v_r}{r}\right) & \dfrac{1}{r}\dfrac{\partial v_z}{\partial \theta}+\dfrac{\partial v_\theta}{\partial z} \\
\dfrac{\partial v_r}{\partial z}+\dfrac{\partial v_z}{\partial r} & \dfrac{1}{r}\dfrac{\partial v_z}{\partial \theta}+\dfrac{\partial v_\theta}{\partial z} & 2\dfrac{\partial v_z}{\partial z}
\end{bmatrix} $ Eq. (03)
Notice that in Eq. (03) $\mathbf{\tau_0}$ is in fact a matrix whose elements are all constants and equal to $\tau_0$. Avoid confusion between the matrix $\mathbf{\tau_0}$ and its elements $\tau_0$. $\tau_0$ is a intrinsic parameter of the fluid featuring the stress that should be surpassed in order to get fluid motion and $\mu_\infty$ is called a plastic viscosity.
- the momentum balance equation.- This is a vectorial partial differential equation with components $(r,\theta,z)$ as follows,
$\rho\left(\frac{\partial v_z}{\partial t}+v_r\frac{\partial v_z}{\partial r}+\frac{v_\theta}{r}\frac{\partial v_z}{\partial \theta}+v_z\frac{\partial v_z}{\partial z}\right)=-\left[\frac{1}{r}\frac{\partial \left(r\tau_{rz}\right)}{\partial r}+\frac{1}{r}\frac{\partial \tau_{\theta z}}{\partial \theta}+\frac{\partial \tau_{zz}}{\partial z}\right]-\frac{\partial p}{\partial z}+\rho g_z $ Eq. (06)
The constitutive and the momentum balance equations form a system of differential equations for the stress $\mathbf{\tau}$ and the fluid velocity $\mathbf{v}$. These equations can be simplified by using the above listed observations to the physics of the problem. First, Eqs. (02-03) are simplified to get,
$\begin{bmatrix}
\tau_{rr} & \tau_{r\theta} & \tau_{rz}\\
\tau_{r\theta} & \tau_{\theta \theta} & \tau_{z\theta}\\
\tau_{rz} & \tau_{z\theta} & \tau_{zz}
\end{bmatrix}
=\begin{bmatrix}
\tau_0 & \tau_0 & \tau_0 \\
\tau_0 & \tau_0 & \tau_0 \\
\tau_0 & \tau_0 & \tau_0
\end{bmatrix}
-\mu_\infty
\begin{bmatrix}
0 & 0 & \dfrac{\partial v_z^{>}}{\partial r} \\
0 & 0 & 0 \\
\dfrac{\partial v_z^{>}}{\partial r} & 0 & 0
\end{bmatrix} $ Eq. (08)
$\tau_{rz}=\tau_0-\mu_\infty \frac{d v_z^{>}\left(r\right)}{d r}$ Eq. (09)
$0=-\frac{1}{r}\frac{d \left(r\tau_{rz}\right)}{d r}-\frac{d p\left(z\right)}{d z}+\rho g_z$ Eq. (10)
$0=-\frac{1}{r}\frac{d}{d r}\left[ r\tau_0-r\mu_\infty\frac{d v_z^{>}}{d r}\right]-\frac{d p\left(z\right)}{d z}-\rho g\cos \chi$ Eq. (11)
\tau_{rr} & \tau_{r\theta} & \tau_{rz}\\
\tau_{r\theta} & \tau_{\theta \theta} & \tau_{z\theta}\\
\tau_{rz} & \tau_{z\theta} & \tau_{zz}
\end{bmatrix}
=\begin{bmatrix}
\tau_0 & \tau_0 & \tau_0 \\
\tau_0 & \tau_0 & \tau_0 \\
\tau_0 & \tau_0 & \tau_0
\end{bmatrix}
-\mu_\infty
\begin{bmatrix}
0 & 0 & \dfrac{\partial v_z^{>}}{\partial r} \\
0 & 0 & 0 \\
\dfrac{\partial v_z^{>}}{\partial r} & 0 & 0
\end{bmatrix} $ Eq. (08)
where the superscripts $<,>$ were introduced to indicate the velocity corresponding to the plug flow from that associated to the parabolic contribution, respectively. From Eq. (08) it can be seen that none of the elements in the stress matrix are zero. However, most elements of the stress matrix do not cause fluid motion but $\tau_{rz}$. A physical interpretation of what is being told in Eq. (08) when $\mathbf{\tau}\geq\mathbf{\tau_0}$ is as follows: due to the orientation of the stress applied on the fluid, the motions or deformations only appear in one direction (as pointed out in the list above). In other words, part of the applied stress is not translated into material deformations but is balanced by the yield stress $\mathbf{\tau_0}$.
Since matrix Eq. (08) gathers a set of equations, the one linking the shear stress and the rate of strain can be taken apart. The resulting equation for the shear stress is,
At this stage, Eqs. (07-09) correspond to the constitutive equation representing the Bingham plastic flowing through the inclined pipe shown in Fig. 1. Since $v_z^{>}$, in Eq. (09), depends only on $r$ partial derivatives notation is not needed anymore. The simplification process to obtain Eqs. (07-09) from Eqs. (02-03) is repeated to simplify the components of the momentum balance Eqs. (04-06). The terms in the radial and angular components, Eqs. (04-05), are zero. From the axial component Eq. (06) only the following remains,
In Eq. (10) $g_z$ represents the $z$ component of the acceleration due to gravity vector $\mathbf{g}$. As mentioned early the fluid travels through an inclined pipe whose axis is parallel to $g_z$ so that $g_z=-g\cos \chi$, as schematically shown in Fig. 3. On the other hand, Eqs. (09-10) can be combined to obtain,
 |
| Fig. 3 Detail of the projection of the acceleration due to gravity on the pipe axis. |
$0=-\frac{1}{r}\frac{d}{d r}\left[ r\tau_0-r\mu_\infty\frac{d v_z^{>}}{d r}\right]-\frac{d p\left(z\right)}{d z}-\rho g\cos \chi$ Eq. (11)
From hereon the velocity profile including the plug flow and parabolic contributions can be calculated by using Eqs. (07,11). First, integration of Eq. (07) is straight so that,
which is in agreement with the sketch in Fig. 2 which shows that $v_z^{<}$ is a straight line perpendicular to the pipe axis. For the parabolic section Eq. (11) can be easily integrated twice to find $v_z^{>}$ if variables are separated. The result is,
$v_z^{>}\left(r\right) = \frac{r^2}{4\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] - \frac{C_2}{\mu_\infty} \ln \lvert r \rvert + \frac{\tau_0 r}{\mu_\infty} + C_3$ Eq. (13)
The constants of integration $C_1$, $C_2$ and $C_3$ in Eqs. (12-13) are unexpected since only the two following common boundary conditions have been considered so far,
$v_z^{>}=0,\quad\text{at}\quad r= R$ Eq. (14)
Evaluation of Eq. (14) are enough to find $C_2$ and $C_3$, but $C_1$ would remain unknown. In other words, a third boundary condition is required. The fact that both, the plug and the parabolic velocities, should coincide at some value of radius can be used to set the missing condition as,
where $r_0$ is unknown but can be found by considering yet another condition based on Eq. (15). Equation (15) establishes that both profiles are the same but smoothness of the change from $v_z^{>}$ to $v_z^{<}$, or vice versa, should be assured for physical reasons. A sudden change from parabolic velocity $\left(v_z^{>}\right)$ to plug velocity $\left(v_z^{<}\right)$, or vice versa, may not correspond to a real life Bingham plastic fluid. Finally, smoothness can be achieved if the derivatives of $v_z^{>}$ and $v_z^{<}$ are the same at $r=r_0$,
Next, $C_2$ and $C_3$ in Eq. (13) are now determined to be,
$C_3=-\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right]\frac{R^2}{4\mu_\infty} - \frac{\tau_0}{\mu_\infty}R$ Eq. (17)
so that $v_z^{>}$ can now be rewritten as follows,
$v_z^{>} \left(r\right) = \frac{R^2}{4\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] \left[\left(\frac{r}{R}\right)^2-1\right] + \frac{\tau_0 R}{\mu_\infty}\left( \frac{r}{R}-1 \right)$ Eq. (18)
The remaining constant $C_1$ is calculated by substituting Eq. (12,18}) into condition Eq. (15). Thus,
Evaluation of the smoothness condition Eq. (16) gives the following relationship for $r_0$,
Equation (20) is a very important result since it defines the position at which the plug and the parabolic velocity are the same. However, despite $r_0$ is a constant it is not numerically defined but in terms of physical parameters or flow conditions. Notice that $r_0$ depends on the pressure gradient, which is implicit in the explanation given at the beginning of this solution, and on $\tau_0$.Also, since $\tau_0$ is a physical parameter of the fluid it would be incorrect to determined it from Eq. (20) instead of using experimental measurements.
A more useful perspective of $r_0$ can be obtained if the pressure terms in Eq. (20) are expressed as shear stress effects. If $\tau_{rz}$ given in Eq. (09) is combined with $v_z^{>}$ in Eq. (18) and then evaluated at the pipe wall $r=R$ the following constant relationship is obtained,
Further combination of Eqs. (20-21) produces a new form of $r_0$,
Equation (22) basically establishes an equivalence between the ratios of radii and shear stresses. However, based on the previous comments and findings, it seems more convenient to write the velocity profile in terms of $\tau_0$ and $\tau_w$ than in terms of $r_0$ and $R$ because speaking of the shear stress is more meaningful and practical. If you know the experimental value of $\tau_0$, an estimation of the velocity or flow rate can be easily made while $r_0$ offers little advantages. Then, $v_z^{<}$ in Eq. (19) is rewritten as,
$v_z^{<} \left(r\right)= \frac{R^2}{4\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] \left[\left(\frac{\tau_0}{\tau_w}\right)^2-1\right] - \left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right]\frac{\tau_0 R^2}{2\tau_w\mu_\infty}\left( \frac{\tau_0}{\tau_w}-1 \right)$
and collecting common factors,
$v_z^{<} \left(r\right)= -\frac{R^2}{4\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] \left[1- \frac{2\tau_0}{\tau_w}+\left(\frac{\tau_0}{\tau_w}\right)^2\right]$ Eq. (23)
A similar process can be performed on Eq. (18) to introduce $\tau_w$ and factorize the pressure terms. Then,
$v_z^{>} \left(r\right) = \frac{R^2}{4\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] \left\lbrace \left[\left(\frac{r}{R}\right)^2-1\right] - \frac{2\tau_0 }{\tau_w }\left( \frac{r}{R}-1 \right) \right\rbrace$ Eq. (24)
Equations (23-24) represent the velocity profile for a Bingham plastic fluid traveling through an inclined pipe. However, $v_z^{<}$ and $v_z^{>}$ are unfinished because the pressure gradient remains implicit and it can not determined at this point due to the piece-wise nature of the velocity profile. Despite this temporal inconvenient the flow rate can still be calculated.
The volumetric flow rate can be calculated by integrating $v_z^{<}$ and $v_z^{>}$ over the cross-section area of the pipe as follows,
$Q=2\pi\int_{0}^{r_0} v_z^{<} rdr + 2\pi\int_{r_0}^{R} v_z^{>} rdr, \quad\text{introduce change of variable}\quad w=\frac{r}{R}$
$Q=-\left. \frac{\pi R^4}{4\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] \left[1- \frac{2\tau_0}{\tau_w}+\left(\frac{\tau_0}{\tau_w}\right)^2\right] \left(\frac{r}{R}\right)^2\right|_{0}^{r_0}$
$+ \left. \frac{\pi R^4}{2\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] \left\lbrace \left[\frac{1}{4}\left(\frac{r}{R}\right)^4-\frac{1}{2}\left(\frac{r}{R}\right)^2\right] - \frac{2\tau_0 }{\tau_w }\left[ \frac{1}{3}\left(\frac{r}{R}\right)^3-\frac{1}{2}\left(\frac{r}{R}\right)^2 \right] \right\rbrace\right|_{r0}^{R}$
after evaluation of the limits some simplifications are made as it follows,
$+ \frac{\pi R^4}{2\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] \left\lbrace -\frac{1}{4} + \frac{\tau_0 }{3\tau_w } - \left[\frac{1}{4}\left(\frac{r_0}{R}\right)^4-\frac{1}{2}\left(\frac{r_0}{R}\right)^2\right] + \frac{2\tau_0 }{\tau_w }\left[ \frac{1}{3}\left(\frac{r_0}{R}\right)^3-\frac{1}{2}\left(\frac{r_0}{R}\right)^2 \right] \right\rbrace$
$Q=\frac{\pi R^4}{4\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] \left\lbrace -\left[1- \frac{2\tau_0}{\tau_w}+\left(\frac{\tau_0}{\tau_w}\right)^2\right] \left(\frac{r_0}{R}\right)^2 -\frac{1}{2} + \frac{2\tau_0 }{3\tau_w } - \left[\frac{1}{2}\left(\frac{r_0}{R}\right)^4-\left(\frac{r_0}{R}\right)^2\right] \right.$
$ + \frac{\tau_0 }{\tau_w }\left[ \frac{4}{3}\left(\frac{r_0}{R}\right)^3-2\left(\frac{r_0}{R}\right)^2 \right] \Bigg\} \quad \text{make use of Eq. (22)}$
$Q=-\frac{\pi R^4}{8\mu_\infty}\left[\frac{d p\left(z\right)}{d z}+\rho g\cos \chi\right] \left\lbrace 1 - \frac{4}{3}\frac{\tau_0 }{\tau_w }+\frac{1}{3}\left(\frac{\tau_0}{\tau_w}\right)^4 \right\rbrace$ Eq. (25)
As before the flow rate $Q$ in Eq. (25) is still implicit because the pressure gradient has not been determined yet. One should notice that Eq. (25) is a first order ordinary differential equation for $p$ and it can be easily solved by separation of variables. The result is,
As usual, $C_4$ can be found with help of the corresponding boundary conditions for $p(z)$,
$p=P_0,\quad\text{at}\quad z=0$
$p=P_L,\quad\text{at}\quad z=L$ Eq. (27)
Despite the two available equations for $C_4$ obtained after the evaluation of Eqs. (27) advantage can be taken from this situation. Combination of the two evaluated expressions provides,
so that,
$\frac{d p\left(z\right)}{d z}=-\frac{\left(P_0-P_L\right)}{L}$ Eq. (29)
Further substitution of Eq. (29) into Eq. (25) gives,
which is an explicit expression for the flow rate $Q$. Some comments are now made on Eq. (30),
- if $\tau_0$ is set to zero the flow rate equation for a newtonian fluid is recovered,
- notice that for any positive value of $\tau_0$ the quantity in braces is always positive too,
- notice that the pressure terms in brackets should always be positive since the fluid is moving upwards.
On the other hand, $v_z^{<}$ and $v_z^{>}$ should be expressed in terms of $P_0$ and $P_L$ as well. This is the end of the problem solution.
Watch the YouTube video for a more dynamic explanation of the problem solution
Any question? Write in the comments and I shall try to help.
Other stuff of interest
- The Poiseuille flow of a power law fluid through an inclined pipe
- The flow of a newonian fluid through an inclined pipe
- The upward flow of a newtonian fluid through the annular region of two vertical pipes
- Understanding the mechanics of the Bernoulli equation
- Practical situation. A polymeric fluid flowing through a pipe
- Practical situation. A polyisoprene solution flowing through a pipe
=========
Ildebrando.
No comments:
Post a Comment