Next: Dimensional1    Previous: Differential2

Chapter Dif (continue)

8.5 Derivations of the Momentum Equation

Fig. 8.8 The shear stress at different surfaces. All shear stress shown in surface $x$ and $x+dx$.

Previously it was shown that equation \eqref{mom:eq:gov} is equivalent to Newton second law for fluids. Equation \eqref{mom:eq:gov} is also applicable for the small infinitesimal cubic. One direction of the vector equation will be derived for $x$ Cartesian coordinate (see Figure 8.8). Later Newton second law will be used and generalized. For surface forces that acting on the cubic are surface forces, gravitation forces (body forces), and internal forces. The body force that acting on infinitesimal cubic in $x$ direction is \begin{align} \label{dif:eq:momentumG} \widehat{i} \cdot \pmb{f}_{B} = {\pmb{f}_{B}}_x \, dx\,dy\,dz \end{align} The dot product yields a force in the directing of $x$. The surface forces in $x$ direction on the $x$ surface on are \begin{align} \label{dif:eq:1momentumXX} f_{xx} = \left. \tau_{xx} \right|_{x+dx} \times \overbrace{dy\,dz}^{dA_x} - \left. \tau_{xx} \right|_{x} \times \overbrace{dy\,dz}^{dA_x} \end{align} The surface forces in $x$ direction on the $y$ surface on are \begin{align} \label{dif:eq:1momentumXY} f_{xy} = \left. \tau_{yx} \right|_{y+dy} \times \overbrace{dx\,dz}^{dA_y} - \left. \tau_{yx} \right|_{y} \times \overbrace{dx\,dz}^{dA_y} \end{align} The same can be written for the $z$ direction. The shear stresses can be expanded into Taylor series as \begin{align} \label{dif:eq:seriesTaylor} \left. \tau_{ix} \right|_{i+di} = \tau_{ix} + \left. \dfrac{\partial \left( \tau_{ix} \right) }{\partial i} \right|_{i} di+ \cdots \end{align} where $i$ in this case is $x$, $y$, or $z$. Hence, the total net surface force results from the shear stress in the $x$ direction is \begin{align} \label{dif:eq:1momentumT} f_x = \left(\dfrac{\partial \tau_{xx} }{\partial x\dfrac{}{}} + \dfrac{\partial \tau_{yx} }{\partial y} + \dfrac{\partial \tau_{zx} }{\partial z} \right) \,dx\, dy\,dz \end{align} after rearrangement equations such as \eqref{dif:eq:1momentumXX} and \eqref{dif:eq:1momentumXY} transformed into \begin{multline} \label{dif:eq:momentum-x} \overbrace{\dfrac{DU_x}{Dt}\,\rho\,\cancel{dx}\,\cancel{dy}\,\cancel{dz}} ^{\text{ internal forces}} = \overbrace{ \left( \dfrac{\partial \tau_{xx} }{\partial x \dfrac{}{}} + \dfrac{\partial \tau_{yx} }{\partial y} + \dfrac{\partial \tau_{zx} }{\partial z} \right) \,\cancel{dx}\, \cancel{dy}\,\cancel{dz} }^{\text{surface forces}} \\ + \overbrace{{f_{G}}_x \,\rho\,\cancel{dx}\, \cancel{dy}\,\cancel{dz}}^{\text{body forces}} \end{multline} equivalent equation \eqref{dif:eq:momentum-x} for $y$ coordinate is \begin{align} \label{dif:eq:momentum-y} \rho\,\dfrac{DU_y}{Dt} = \dfrac{\partial \tau_{xy} }{\partial x} + \dfrac{\partial \tau_{yy} }{\partial y} + \dfrac{\partial \tau_{zy} }{\partial z} + \rho\,{f_{G}}_y \end{align} The same can be obtained for the $z$ component \begin{align} \label{dif:eq:momentum-z} \rho\,\dfrac{DU_z}{Dt} = \dfrac{\partial \tau_{xz} }{\partial x} + \dfrac{\partial \tau_{yz} }{\partial y} + \dfrac{\partial \tau_{zz} }{\partial z} + \rho\,{f_{G}}_z \end{align}

Generally the component momentum equation is as \begin{align} \label{dif:eq:momentum-i} \rho \, \dfrac{DU_i}{Dt} = \dfrac{\partial \tau_{ii} }{\partial i} + \dfrac{\partial \tau_{ji} }{\partial j} + \dfrac{\partial \tau_{ki} }{\partial j} + \rho \, {f_{G}}_i \end{align}

Where $i$ is the balance direction and $j$ and $k$ are two other coordinates. Equation \eqref{dif:eq:momentum-i} can be written in a vector form which combined all three components into one equation. The advantage of the vector from allows the usage of the different coordinates. The vector form is \begin{align} \label{dif:eq:momentum-V} \rho \, \dfrac{D\pmb{U}}{Dt} = \nabla \cdot \boldsymbol{\tau^{(i)}} + \rho \, {\pmb{f_{G}}} \end{align} where here \begin{align} \nonumber \boldsymbol{\tau^{(i)}} = \tau_{ix} \widehat{i} + \tau_{iy} \widehat{j} + \tau_{iz} \widehat{k} \end{align} is part of the shear stress tensor and $i$ can be any of the $x, y,$ or $z$. Or in index (Einstein) notation as \begin{align} \label{dif:eq:momentum-E} \rho \, \dfrac{DU_i}{Dt} = \dfrac{\partial \tau_{ji} }{\partial x_i} + \rho \, {f_{G}}_i \end{align}

Equations \eqref{dif:eq:momentum-x} or \eqref{dif:eq:momentum-y} or \eqref{dif:eq:momentum-z} requires that the stress tensor be defined in term of the velocity/deformation. The relationship between the stress tensor and deformation depends on the classes of materials the stresses acts on. Additionally, the deformation can be viewed as a function of the velocity field. As engineers do in general, the simplest model is assumed which referred as the solid continuum model. In this model the relationship between the (shear) stresses and rate of strains are assumed to be linear. In solid material, the shear stress yields a fix amount of deformation. In contrast, when applying the shear stress in fluids, the result is a continuous deformation. Furthermore, reduction of the shear stress does not return the material to its original state as in solids. The similarity to solids the increase shear stress in fluids yields larger deformations. Thus this solid'' model is a linear relationship with three main assumptions:

Fig. 8.9 Control volume at $t$ and $t + dt$ under continuous angle deformation. Notice the three combinations of the deformation shown by purple color relative to blue color.

1. There is no preference in the orientation (also call isentropic fluid),
2. there is no left over stresses (In other words when the no shear stress'' situation exist the rate of deformation or strain is zero), and
3. a linear relationship exist between the shear stress and the rate of shear strain.
At time $t$, the control volume is at a square shape and at a location as depicted in Figure 8.9 (by the blue color). At time $t +dt$ the control volume undergoes three different changes. The control volume moves to a new location, rotates and changes the shape (the purple color in in Figure 8.9). The translational movement is referred to a movement of body without change of the body and without rotation. The rotation is the second movement that referred to a change in of the relative orientation inside the control volume. The third change is the misconfiguration or control volume (deformation). The deformation of the control volume has several components (see the top of Figure 8.9). The shear stress is related to the change in angle of the control volume lower left corner. The angle between $x$ to the new location of the control volume can be approximate for a small angle as \begin{align} \label{dif:eq:dGammax} \dfrac{d\gamma_x}{dt} = \tan\left( \dfrac{U_y + \dfrac{dU_y}{dx}dx - U_y } {dx} \right) = \tan\left( \dfrac{dU_y}{dx} \right) \cong \dfrac{dU_y}{dx} \end{align} The total angle deformation (two sides $x$ and $y$) is \begin{align} \label{dif:eq:shearAngle} \dfrac{D\gamma_{xy}}{Dt} = \dfrac{dU_y}{dx} + \dfrac{dU_x}{dy} \end{align} In these derivatives, the symmetry $\dfrac{dU_y}{dx} \neq \dfrac{dU_x}{dy}$ was not assumed and or required because rotation of the control volume. However, under isentropic material it is assumed that all the shear stresses contribute equally. For the assumption of a linear fluid . \begin{align} \label{dif:eq:linearFluidXY} \tau_{xy} = \mu \dfrac{D\gamma_{xy}}{Dt} = \mu \left( \dfrac{dU_y}{dx} + \dfrac{dU_x}{dy} \right) \end{align} where, $\mu$ is the normal'' or ordinary'' viscosity coefficient which relates the linear coefficient of proportionality and shear stress. This deformation angle coefficient is assumed to be a property of the fluid. In a similar fashion it can be written to other directions for $x\,z$ as \begin{align} \label{dif:eq:linearFluidXZ} \tau_{xz} = \mu \dfrac{D\gamma_{xz}}{Dt} = \mu \left( \dfrac{dU_z}{dx} + \dfrac{dU_x}{dz} \right) \end{align} and for the directions of $y\,z$ as \begin{align} \label{dif:eq:linearFluidYZ} \tau_{yz} = \mu \dfrac{D\gamma_{yz}}{Dt} = \mu \left( \dfrac{dU_z}{dy\dfrac{}{}} + \dfrac{dU_y}{dz} \right) \end{align} is assumed to be the same regardless of the direction. This assumption is referred as isotropic viscosity. It can be noticed at this stage, the relationship for the two of stress tensor parts was established. The only missing thing, at this stage, is the diagonal component which to be dealt below.

In general equation \eqref{dif:eq:linearFluidXY} can be written as \begin{align} \label{dif:eq:linearFluidIJ} \tau_{ij} = \mu \dfrac{D\gamma_{ij}}{Dt} = \mu \left( \dfrac{dU_j}{di\dfrac{}{}} + \dfrac{dU_i}{dj} \right) \end{align} where $i\neq j$ and $i=x$ or $y$ or $z$.

Normal Stress

The normal stress, $\tau_{ii}$ (where $i$ is either ,$x$, $y$, $z$) appears in the shear matrix diagonal. To find the main (or the diagonal) stress the coordinates are rotate by $45^\circ$. The diagonal lines (line and line in Figure ) in the control volume move to the new locations. In addition, the sides and rotate in unequal amount which make one diagonal line longer and one diagonal line shorter. The normal shear stress relates to the change in the diagonal line length change. This relationship can be obtained by changing the coordinates orientation as depicted by Figure 8.10. The $dx$ is constructed so it equals to $dy$. The forces acting in the direction of $x'$ on the element are combination of several terms. For example, on the $x$'' surface (lower surface) and the $y$'' (left) surface, the shear stresses are acting in this direction. It can be noticed that $d x'$'' surface is $\sqrt{2}$ times larger than $dx$ and $dy$ surfaces. The force balance in the $x'$ is \begin{align} \label{dif:eq:tauxxxpxpT} \overbrace{dy}^{A_x} \tau_{xx} \overbrace{\dfrac{1}{\sqrt{2}} }^{\cos\theta_x} + \overbrace{dx}^{A_y} \tau_{yy} \overbrace{\dfrac{1}{\sqrt{2}} }^{\cos\theta_y} + \overbrace{dx}^{A_y} \tau_{yx} \overbrace{\dfrac{1}{\sqrt{2}} }^{\cos\theta_y} + \overbrace{dy}^{A_x} \tau_{xy} \overbrace{\dfrac{1}{\sqrt{2}} }^{\cos\theta_y} = \overbrace{dx\sqrt{2}}^{A_{x'} } \,\tau_{ x' x' } \end{align} dividing by $dx$ and after some rearrangements utilizing the identity $\tau_{xy}= \tau_{yx}$ results in \begin{align} \label{dif:eq:tauxxxpxp} \dfrac{\tau_{xx} + \tau_{yy}}{2} + \tau_{yx} = \tau_{x'x'} \end{align} Setting the similar analysis in the $y'$ results in \begin{align} \label{dif:eq:tauyyypyp} \dfrac{\tau_{xx} + \tau_{yy}}{2} - \tau_{yx} = \tau_{y'y'} \end{align} Subtracting \eqref{dif:eq:tauyyypyp} from \eqref{dif:eq:tauxxxpxp} results in \begin{align} \label{dif:eq:tauxpxpypyp} 2\,\tau_{yx} = \tau_{x'x'} - \tau_{y'y'} \end{align}

Figure 8.11 Different triangles deformation for the calculations of the normal stress.

or dividing by $2$ equation \eqref{dif:eq:tauxpxpypyp} becomes \begin{align} \label{dif:eq:tauxpxpypypD} \tau_{yx} = \dfrac{1}{2} \left( \tau_{x' x'} - \tau_{y' y'} \right) \end{align} Equation \eqref{dif:eq:tauxpxpypyp} relates the difference between the normal shear stress and the normal shear stresses in $x'$, $y'$ coordinates) and the angular strain rate in the regular ($x, y$ coordinates). The linear deformations in the $x'$ and $y'$ directions which is rotated $45^{\circ}$ relative to the $x$ and $y$ axes can be expressed in both coordinates system. The angular strain rate in the ($x$, $y$) is frame related to the strain rates in the ($x'$, $y'$) frame. Figure 8.11(a) depicts the deformations of the triangular particles between time $t$ and $t+dt$. The small deformations a , b, c, and d in the Figure are related to the incremental linear strains. The rate of strain in the $x$ direction is \begin{align} \label{dif:eq:difEx} d\epsilon_x = \dfrac{c}{dx} \end{align} The rate of the strain in $y$ direction is \begin{align} \label{dif:eq:difEy} d\epsilon_y = \dfrac{a}{dx} \end{align} The total change in the deformation angle is related to $\tan\theta$, in both sides ($d/dx + b/dy$) which in turn is related to combination of the two sides angles. The linear angular deformation in $xy$ direction is \begin{align} \label{dif:eq:difGy} d\gamma_{xy} = \dfrac{b+d}{dx} \end{align} Here, $d\epsilon_x$ is the linear strain (increase in length divided by length) of the particle in the $x$ direction, and $d\epsilon_y$ is its linear strain in the y-direction. The linear strain in the $x'$ direction can be computed by observing Figure 8.11(b). The hypotenuse of the triangle is oriented in the $x'$ direction (again observe Figure 8.11(b)). The original length of the hypotenuse $\sqrt{2} dx$. The change in the hypotenuse length is $\sqrt{\left( c+b\right)^2 + \left( a+d\right)^2}$. It can be approximated that the change is about $45^{\circ}$ because changes are infinitesimally small. Thus, $\cos 45^{\circ}$ or $\sin 45^{\circ}$ times the change contribute as first approximation to change. Hence, the ratio strain in the $x'$ direction is \begin{align} \label{dif:eq:difExp} d\epsilon_{x'} = \dfrac{\sqrt{\left( c+b\right)^2 + \left( a+d\right)^2}}{\sqrt{2}dx} \simeq \dfrac{ \dfrac{\left( c+b\right)}{\sqrt{2}} + \dfrac{\left( c+b\right)}{\sqrt{2}} + \overbrace{f \left( d x'\right)}^{\sim 0} } {\sqrt{2}dx} \end{align} Equation \eqref{dif:eq:difExp} can be interpreted as (using equations \eqref{dif:eq:difEx}, \eqref{dif:eq:difEy}, and \eqref{dif:eq:difGy}) as \begin{align} \label{dif:eq:difExpi} d\epsilon_{x'} = \dfrac{1}{2} \left(\dfrac{a+b+c+d}{dx\dfrac{}{}} \right) = \dfrac{1}{2} \left(d\epsilon_y + d\epsilon_y + d\gamma_{xy} \right) \end{align} In the same fashion, the strain in $y'$ coordinate can be interpreted to be \begin{align} \label{dif:eq:difEypi} d\epsilon_{y'} = \dfrac{1}{2} \left(d\epsilon_y + d\epsilon_y - d\gamma_{xy} \right) \end{align} Notice the negative sign before $d\gamma_{xy}$. Combining equation \eqref{dif:eq:difExpi} with equation \eqref{dif:eq:difEypi} results in \begin{align} \label{dif:eq:defExyComb} d\epsilon_{x'} - d\epsilon_{y' } = d\gamma_{xy} \end{align} Equation \eqref{dif:eq:defExyComb} describing in Lagrangian coordinates a single particle. Changing it to the Eulerian coordinates transforms equation qref{dif:eq:defExyComb} into \begin{align} \label{dif:eq:defExyCombE} \dfrac{D\epsilon_{x'}}{Dt} - \dfrac{D\epsilon_{y'}}{Dt} = \dfrac{D\gamma_{xy}}{Dt} \end{align} From \eqref{dif:eq:linearFluidXY} it can be observed that the right hand side of equation \eqref{dif:eq:defExyCombE} can be replaced by $\tau_{xy}/\mu$ to read \begin{align} \label{dif:eq:defExyCombEIntermindateStep} \dfrac{D\epsilon_{x'}}{Dt} - \dfrac{D\epsilon_{y'}}{Dt} = \dfrac{\tau_{xy}}{\mu} \end{align} From equation \eqref{dif:eq:tauxpxpypyp} $\tau_{xy}$ be substituted and equation \eqref{dif:eq:defExyCombEIntermindateStep} can be continued and replaced as \begin{align} \label{dif:eq:tauExyCombE} \dfrac{D\epsilon_{x'}}{Dt} - \dfrac{D\epsilon_{y'}}{Dt} = \dfrac{1}{2\,\mu}\left( \tau_{x'x'} - \tau_{y'y'} \right) \end{align}

Fig. 8.12 Linear strain of the element purple denotes $t$ and blue is for $t+dt$. Dashed squares denotes the movement without the linear change.

Figure 8.12 depicts the approximate linear deformation of the element. The linear deformation is the difference between the two sides as \begin{align} \label{dif:eq:linearDeformation} \dfrac {D\epsilon_{x'} }{Dt} = \dfrac {\partial U_{x'} } {\partial x' } \end{align} The same way it can written for the $y{\kern -1pt \lower+1pt\hbox{'}}$ coordinate. \begin{align} \label{dif:eq:linearDeformationYP} \dfrac {D\epsilon_{y'} }{Dt} = \dfrac {\partial U_{y'} } {\partial y' } \end{align} Equation \eqref{dif:eq:linearDeformation} can be written in the $y{\kern -1pt \lower+1pt\hbox{'}}$ and is similar by substituting the coordinates. The rate of strain relations can be substituted by the velocity and equations \eqref{dif:eq:linearDeformation} and \eqref{dif:eq:linearDeformationYP} changes into \begin{align} \label{dif:eq:almostMufractionxpyp} \tau_{x'x'} - \tau_{y'y'} = 2 \mu \left( \dfrac {\partial U_{x' } } {\partial x' } - \dfrac {\partial U_{y' } } {\partial y' } \right) \end{align} Similar two equations can be obtained in the other two plans. For example in $y\kern -1pt \lower+1pt\hbox{'}$– $z\kern -1pt \lower+1pt\hbox{'}$ plan one can obtained \begin{align} \label{dif:eq:almostMufractionypzp} \tau_{x'x'} - \tau_{z'z'} = 2 \mu \left( \dfrac {\partial U_{x' } } {\partial x' } - \dfrac {\partial U_{z' } } {\partial z' } \right) \end{align} Adding equations \eqref{dif:eq:almostMufractionxpyp} and \eqref{dif:eq:almostMufractionypzp} results in \begin{align} \label{dif:eq:totalStressDefraction} \overbrace{\left( 3 - 1\right)}^{2} \, \tau_{x'x'} - \tau_{y'y'} - \tau_{z'z'} = \overbrace{\left( 6- 2 \right)}^{4} \mu \dfrac {\partial U_{x' } } {\partial x' } - 2\,\mu \left( \dfrac {\partial U_{y' } } {\partial y'\dfrac{}{} } + \dfrac {\partial U_{z' } } {\partial z' } \right) \end{align} rearranging equation \eqref{dif:eq:totalStressDefraction} transforms it into \begin{align} \label{dif:eq:totalStressDefraction1} 3\, \tau_{x'x'} = \tau_{x'x'} + \tau_{y'y'} + \tau_{z'z'} + {6}\, \mu \dfrac {\partial U_{x'} } {\partial x' } - 2\,\mu \left( \dfrac {\partial U_{x' } }{\partial x'\dfrac{}{} } + \dfrac {\partial U_{y' } }{\partial y' } + \dfrac {\partial U_{z' } }{\partial z' } \right) \end{align} Dividing the results by 3 so that one can obtained the following \begin{align} \label{dif:eq:totalStressDefraction2} \tau_{x'x'} = \overbrace{ \dfrac{\tau_{x'x'} + \tau_{y'y'} + \tau_{z'z'} }{3} }^{\text{mechanical'' pressure} } + {2}\, \mu \dfrac {\partial U_{x' } }{\partial x' } - \dfrac{2}{3}\,\mu \left( \dfrac {\partial U_{x' } }{\partial x'\dfrac{}{} } + \dfrac {\partial U_{y' } }{\partial y' } + \dfrac {\partial U_{z' } }{\partial z' } \right) \end{align} The mechanical'' pressure, $P_m$, is defined as the (negative) average value of pressure in directions of $x'$–$y'$–$z'$. This pressure is a true scalar value of the flow field since the propriety is averaged invariant to the coordinate transformation. In situations where the main diagonal terms of the stress tensor are not the same in all directions (in some viscous flows) this property can be served as a measure of the local normal stress. The mechanical pressure can be defined as averaging of the normal stress acting on a infinitesimal sphere. It can be shown that this two definitions are identical'' in the limits . With this definition and noticing that the coordinate system $x'$–$y'$ has no special significance and hence equation \eqref{dif:eq:totalStressDefraction2} must be valid in any coordinate system thus equation \eqref{dif:eq:totalStressDefraction2} can be written as \begin{align} \label{dif:eq:tStressDeft} \tau_{xx} = - P_m + 2 \, \mu \dfrac{\partial U_x}{\partial x} + \dfrac{2}{3} \mu \, \nabla \cdot \pmb{U} \end{align} Again where $P_m$ is the mechanical pressure and is defined as

Mechanical Pressure

\begin{align} \label{dif:eq:mechanicalP} P_m = - \dfrac{\tau_{xx} + \tau_{yy} + \tau_{zz} }{3} \end{align}
It can be observed that the non main (diagonal) terms of the stress tensor are represented by an equation like \eqref{dif:eq:linearFluidIJ}. Commonality engineers like to combined the two difference expressions into one as \begin{align} \label{dif:eq:combinedStress3} \tau_{xy} = - \left( P_m + \dfrac{2}{3\dfrac{}{}}\mu \nabla \cdot \pmb{U} \right) \overbrace{\delta_{xy}}^{=0} + \mu \left( \dfrac{\partial U_x}{\partial y \dfrac{}{}} + \dfrac{\partial U_y}{\partial x } \right) \end{align} or \begin{align} \label{dif:eq:combinedStress2} \tau_{xx} = - \left( P_m + \dfrac{2}{3\dfrac{}{}}\mu \nabla \cdot \pmb{U} \right) \overbrace{\delta_{xy}}^{=1} + \mu \left( \dfrac{\partial U_x}{\partial x\dfrac{}{} } + \dfrac{\partial U_y}{\partial y } \right) \end{align}

or index notation \begin{align} \label{dif:eq:combinedStress1} \tau_{ij} = - \left( P_m + \dfrac{2}{3\dfrac{}{}}\mu \nabla \cdot \pmb{U} \right) \delta_{ij} + \mu \left( \dfrac{\partial U_i}{\partial x_j\dfrac{}{}} + \dfrac{\partial U_j}{\partial x_i } \right) \end{align}

where $\delta_{ij}$ is the Kronecker delta what is $\delta_{ij}= 1$ when $i=j$ and $\delta_{ij}= 0$ otherwise. While this expression has the advantage of compact writing, it does not add any additional information. This expression suggests a new definition of the thermodynamical pressure is

Thermodynamic Pressure

\begin{align} \label{dif:eq:thermoP} P = P_m + \dfrac{2}{3}\mu \nabla \cdot \pmb{U} \end{align}
}

Summary of The Stress Tensor

The above derivations were provided as a long mathematical explanation . To reduced one unknown (the shear stress) equation \eqref{dif:eq:momentum-x} the relationship between the stress tensor and the velocity were to be established. First, connection between $\tau_{xy}$ and the deformation was built. Then the association between normal stress and perpendicular stress was constructed. Using the coordinates transformation, this association was established. The linkage between the stress in the rotated coordinates to the deformation was established.

Second Viscosity Coefficient

The coefficient ${2}/{3}\mu$ is experimental and relates to viscosity. However, if the derivations before were to include additional terms, an additional correction will be needed. This correction results in \begin{align} \label{dif:eq:thermoPLemda} P = P_m + \lambda \nabla \cdot \pmb{U} \end{align} The value of $\lambda$ is obtained experimentally. This coefficient is referred in the literature by several terms such as the expansion viscosity'' second coefficient of viscosity'' and bulk viscosity.'' Here the term bulk viscosity will be adapted. The dimension of the bulk viscosity, $\lambda$, is similar to the viscosity $\mu$. According to second law of thermodynamic derivations (not shown here and are under construction) demonstrate that $\lambda$ must be positive. The thermodynamic pressure always tends to follow the mechanical pressure during a change. The expansion rate of change and the fluid molecular structure through $\lambda$ control the difference. Equation \eqref{dif:eq:thermoPLemda} can be written in terms of the thermodynamic pressure $P$, as \begin{align} \label{did:eq:combinedStress} \tau_{ij} = - \left[ P + \left( \dfrac{2}{3\dfrac{}{}}\mu - \lambda \right) \nabla \cdot \pmb{U} \right] \delta_{ij} + \mu \left( \dfrac{\partial U_i}{\partial x_j\dfrac{}{}} + \dfrac{\partial U_j}{\partial x_i } \right) \end{align} The significance of the difference between the thermodynamic pressure and the mechanical pressure associated with fluid dilation which connected by $\nabla \cdot \pmb{U}$. The physical meaning of $\nabla \cdot \pmb{U}$ represents the relative volume rate of change. For simple gas (dilute monatomic gases) it can be shown that $\lambda$ vanishes. In material such as water, $\lambda$ is large (3 times $\mu$) but the net effect is small because in that cases $\nabla \cdot \pmb{U}\longrightarrow 0$. For complex liquids this coefficient, $\lambda$, can be over 100 times larger than $\mu$. Clearly for incompressible flow, this coefficient or the whole effect is vanished . In most cases, the total effect of the dilation on the flow is very small. Only in micro fluids and small and molecular scale such as in shock waves this effect has some significance. In fact this effect is so insignificant that there is difficulty in to construct experiments so this effect can be measured. Thus, neglecting this effect results in \begin{align} \label{did:eq:combinedStressNoLambda} \tau_{ij} = - P \delta_{ij} + \mu \left( \dfrac{\partial U_i}{\partial x_j\dfrac{}{}} + \dfrac{\partial U_j}{\partial x_i } \right) \end{align} To explain equation \eqref{did:eq:combinedStressNoLambda}, it can be written for specific coordinates. For example, for the $\tau_{xx}$ it can be written that \begin{align} \label{dif:eq:combinedStressNoLambda–xx} \tau_{xx} = - P + 2 \dfrac{\partial U_x}{\partial x } \end{align} and the $y$ coordinate the equation is \begin{align} \label{dif:eq:combinedStressNoLambda–yy} \tau_{yy} = - P + 2 \dfrac{\partial U_y}{\partial y } \end{align} however the mix stress, $\tau_{xy}$, is \begin{align} \label{dif:eq:combinedStressNoLambda–xy} \tau_{xy} = \tau_{yx} = \left( \dfrac{\partial U_y}{\partial x\dfrac{}{} } + \dfrac{\partial U_x}{\partial y } \right) \end{align} For the total effect, substitute equation \eqref{did:eq:combinedStress} into equation \eqref{dif:eq:momentum-x} which results in \begin{align} \label{dif:eq:nsGx} \begin{array}{rrcll} \rho\left( \dfrac{D U_x}{Dt\dfrac{}{}} \right) = -\dfrac{\partial }{\partial x} \left( P+ \left(\dfrac{2}{3\dfrac{}{}}\mu - \lambda \right) \nabla\cdot \pmb{U} \right) + \mu\,\left( \dfrac{\partial^2 U_{x} }{\partial x^2 \dfrac{}{}} + \dfrac{\partial^2 U_{x} }{\partial y^2} + \dfrac{\partial^2 U_{x} }{\partial z^2} \right) & + {\pmb{f}_{B}}_{x} \end{array} \end{align} or in a vector form as

N-S in stationary Coordinates

\begin{align} \label{dif:eq:nsGv} \rho \, \dfrac{D \pmb{U}}{Dt} = - \nabla P + \left(\dfrac{1}{3\dfrac{}{}}\mu + \lambda \right) \nabla \, \left( \nabla \cdot \pmb{U} \right) + \mu\,\nabla^{2} \pmb{U} + {\pmb{f}_{B}} \end{align}
For in index form as \begin{align} \label{dif:eq:nsG} \rho \, \dfrac{D\, U_i}{Dt} = - \dfrac{\partial }{\partial x_i} \left( P+ \left(\dfrac{2}{3\dfrac{}{}}\mu - \lambda \right) \nabla\cdot \pmb{U} \right) + \dfrac { \partial }{\partial x_j} \left( \mu \left( \dfrac{\partial U_i }{\partial x_j\dfrac{}{}} + \dfrac{\partial U_j} {\partial x_i} \right) \right) + {\pmb{f}_{B}}_i \end{align} For incompressible flow the term $\nabla\cdot\pmb{U}$ vanishes, thus equation \eqref{dif:eq:nsGv} is reduced to

Momentum for Incompressible Flow

\begin{align} \label{dif:eq:nsGvIncompressibleFlow1} \rho \, \dfrac{D \pmb{U}}{Dt} = - \nabla P + \mu\,\nabla^{2} \pmb{U} + {\pmb{f}_{B}} \end{align}
or in the index notation it is written \begin{align} \label{dif:eq:nsGvIncompressibleFlow} \rho \, \dfrac{D\, U_i}{Dt} = - \dfrac{\partial P }{\partial x_i} + \mu \, \dfrac { \partial^2 \pmb{U} }{\partial x_i \partial x_j} + {\pmb{f}_{B}}_i \end{align} The momentum equation in Cartesian coordinate can be written explicitly for $x$ coordinate as \begin{multline} \label{dif:eq:momEqx} \rho\, \left(\dfrac{\partial U_x}{\partial t\dfrac{}{}} + \right. U_x\dfrac{\partial U_x}{\partial x} + U_y \dfrac{\partial U_x}{\partial y} + \left. U_z \dfrac{\partial U_x}{\partial z\dfrac{}{}}\right) = \\ -\dfrac{\partial P}{\partial x} + \mu \left(\dfrac{\partial^2 U_x}{\partial x^2\dfrac{}{}} + \dfrac{\partial^2 U_x}{\partial y^2} + \dfrac{\partial^2 U_x}{\partial z^2}\right) + \rho\, g_x \end{multline} Where $g_x$ is the the body force in the $x$ direction ($\widehat{i}\cdot\pmb{g}$). In the $y$ coordinate the momentum equation is \begin{multline} \label{dif:eq:momEqy} \rho\,\left(\dfrac{\partial U_y}{\partial t\dfrac{}{}} + \right. U_x \dfrac{\partial U_y}{\partial x} + U_y \dfrac{\partial U_y}{\partial y}+ \left. U_z \dfrac{\partial U_y}{\partial z\dfrac{}{}}\right) = \\ -\dfrac{\partial P}{\partial y} + \mu \left(\dfrac{\partial^2 v}{\partial x^2 \dfrac{}{} } + \dfrac{\partial^2 v}{\partial y^2} + \dfrac{\partial^2 v}{\partial z^2}\right) + \rho\, g_y \end{multline} in $z$ coordinate the momentum equation is \begin{multline} \label{dif:eq:momEqz} \rho\, \left(\dfrac{\partial U_z}{\partial t\dfrac{}{}} + \right. U_x \dfrac{\partial U_z}{\partial x} + U_y \dfrac{\partial U_z}{\partial y}+ \left. U_z \dfrac{\partial U_z}{\partial z\dfrac{}{} }\right) = \\ -\dfrac{\partial P}{\partial z} + \mu \left(\dfrac{\partial^2 U_z}{\partial x^2} + \dfrac{\partial^2 U_z}{\partial y^2} + \dfrac{\partial^2 U_z}{\partial z^2\dfrac{}{}}\right) + \rho\, g_z \end{multline}

8.6 Boundary Conditions and Driving Forces

8.6.1 Boundary Conditions Categories

The governing equations that were developed earlier requires some boundary conditions and initial conditions. These conditions described physical situations that are believed or should exist or approximated. These conditions can be categorized by the velocity, pressure, or in more general terms as the shear stress conditions (mostly at the interface). For this discussion, the shear tensor will be separated into two categories, pressure (at the interface direction) and shear stress (perpendicular to the area). A common velocity condition is that the liquid has the same value as the solid interface velocity. In the literature, this condition is referred as the no slip'' condition. The solid surface is rough thus the liquid participles (or molecules) are slowed to be at the solid surface velocity. This boundary condition was experimentally observed under many conditions yet it is not universal true. The slip condition (as oppose to no slip'' condition) exist in situations where the scale is very small and the velocity is relatively very small. The slip condition is dealing with a difference in the velocity between the solid (or other material) and the fluid media. The difference between the small scale and the large scale is that the slip can be neglected in the large scale while the slip cannot be neglected in the small scale. In another view, the difference in the velocities vanishes as the scale increases.

Fig. 8.13 1–Dimensional free surface describing $\widehat{\pmb{n}}$ and $\widehat{\pmb{t}}$.

Another condition which affects whether the slip condition exist is how rapidly of the velocity change. The slip condition cannot be ignored in some regions, when the flow is with a strong velocity fluctuations. Mathematically the no slip'' condition is written as \begin{align} \label{dif:eq:noSlipCondV} \widehat{\mathbf{t}} \cdot \left( \pmb{U}_{fluid} - \pmb{U}_{boundary} \right) = 0 \end{align} where $\widehat{\mathbf{n}}$ is referred to the area direction (perpendicular to the area see Figure \ref{dif:fig:2DfreeSurface}). While this condition \eqref{dif:eq:noSlipCondV} is given in a vector form, it is more common to write this condition as a given velocity at a certain point such as \begin{align} \label{dif:eq:noSlipCond} U(\ell) = U_{\ell} \end{align} Note, the no slip'' condition is applicable to the ideal fluid (inviscid flows'') because this kind of flow normally deals with large scales. The 'slip' condition is written in similar fashion to equation \eqref{dif:eq:noSlipCondV} as \begin{align} \label{dif:eq:slipCondV} \widehat{\mathbf{t}} \cdot \left( \pmb{U}_{fluid} - \pmb{U}_{boundary} \right) = f(Q,scale, etc) \end{align} As oppose to a given velocity at particular point, a requirement on the acceleration (velocity) can be given in unknown position. The condition \eqref{dif:eq:noSlipCondV} can be mathematically represented in another way for free surface conditions. To make sure that all the material is accounted for in the control volume (does not cross the free surface), the relative perpendicular velocity at the interface must be zero. The location of the (free) moving boundary can be given as $f(\widehat{\pmb{r}},t) =0$ as the equation which describes the bounding surface. The perpendicular relative velocity at the surface must be zero and therefore \begin{align} \label{dif:eq:perpendicularUSurface} \dfrac{Df}{Dt} = 0 \quad \mbox{ on the surface } f (\widehat{\pmb{r}},t) = 0 \end{align} This condition is called the kinematic boundary condition. For example, the free surface in the two dimensional case is represented as $f(t,x,y)$. The condition becomes as \begin{align} \label{dif:eq:2perpendicularUSurface} 0 = \dfrac{\partial f}{\partial t} + U_x\, \dfrac{\partial f}{\partial x} + U_y\, \dfrac{\partial f}{\partial y} \end{align} The solution of this condition, sometime, is extremely hard to handle because the location isn't given but the derivative given on unknown location. In this book, this condition will not be discussed (at least not plane to be written). The free surface is a special case of moving surfaces where the surface between two distinct fluids. In reality the interface between these two fluids is not a sharp transition but only approximation (see for the surface theory). There are situations where the transition should be analyzed as a continuous transition between two phases. In other cases, the transition is idealized an almost jump (a few molecules thickness). Furthermore, there are situations where the fluid (above one of the sides) should be considered as weightless material. In these cases the assumptions are that the transition occurs in a sharp line, and the density has a jump while the shear stress are continuous (in some cases continuously approach zero value). While a jump in density does not break any physical laws (at least those present in the solution), the jump in a shear stress (without a jump in density) does break a physical law. A jump in the shear stress creates infinite force on the adjoin thin layer. Off course, this condition cannot be tolerated since infinite velocity (acceleration) is impossible. The jump in shear stress can appear when the density has a jump in density. The jump in the density (between the two fluids) creates a surface tension which offset the jump in the shear stress. This condition is expressed mathematically by equating the shear stress difference to the forces results due to the surface tension. The shear stress difference is \begin{align} \label{dif:eq:JumpStress} \Delta \boldsymbol{\tau}^{(n)} = 0 = \Delta {\boldsymbol{\tau}^{(n)}}_{\text{ upper surface}} - \Delta {\boldsymbol{\tau}^{(n)}}_{\text{ lower surface}} \end{align} where the index $(n)$ indicate that shear stress are normal (in the surface area). If the surface is straight there is no jump in the shear stress. The condition with curved surface are out the scope of this book yet mathematically the condition is given as without explanation as \begin{align} \label{dif:eq:bcCurvedSurface} \widehat{\pmb{n}}\cdot {\boldsymbol{\tau}^{(n)}} = \sigma \left( \dfrac{1\dfrac{}{}}{R_1} + \dfrac{1}{R_2}\right) \ \widehat{\pmb{t}}\cdot {\boldsymbol{\tau}^{(t)}} = - \widehat{\pmb{t}}\cdot \nabla\sigma \end{align} where $\widehat{\pmb{n}}$ is the unit normal and $\widehat{\pmb{t}}$ is a unit tangent to the surface (notice that direction pointed out of the center'' see Figure \ref{dif:fig:2DfreeSurface}) and $R_1$ and $R_2$ are principal radii. One of results of the free surface condition (or in general, the moving surface condition) is that integration constant is unknown). In same instances, this constant is determined from the volume conservation. In index notation equation \eqref{dif:eq:bcCurvedSurface} is written \begin{align} \label{dif:eq:bcCurvedSurfaceIndex} \tau_{ij}^{(1)}\,n_j + \sigma \, n_i \left( \dfrac{1}{R_1} + \dfrac{1}{R_2}\right) = \tau_{ij}^{(2)}\,n_j \end{align} where $1$ is the upper surface and $2$ is the lower surface. For example in one dimensional \begin{align} \label{dif:eq:nAndT} \begin{array}{rl} \widehat{\pmb{n}} = & \dfrac{\left(-f^\prime{}(x), 1 \right) }{\sqrt{1+ \left(f^\prime{}(x)\right)^2 }} \ \widehat{\pmb{t}} = & \dfrac{\left(1,f^\prime{}(x) \right) }{\sqrt{1+ \left(f^\prime{}(x)\right)^2 }} \end{array} \end{align} the unit vector is given as two vectors in $x$ and $y$ and the radius is given by equation \eqref{intro:eq:radius}. The equation is given by \begin{align} \label{dif:eq:1DfreeSurface} \dfrac{\partial f}{\partial t} + U_x \dfrac{\partial f}{\partial x} = U_y \end{align}

The Pressure Condition

The second condition that commonality prescribed at the interface is the static pressure at a specific location. The static pressure is measured perpendicular to the flow direction. The last condition is similar to the pressure condition of prescribed shear stress or a relationship to it. In this category include the boundary conditions with issues of surface tension which were discussed earlier. It can be noticed that the boundary conditions that involve the surface tension are of the kind where the condition is given on boundary but no at a specific location.

Gravity as Driving Force

The body forces, in general and gravity in a particular, are the condition that given on the flow beside the velocity, shear stress (including the surface tension) and the pressure. The gravity is a common body force which is considered in many fluid mechanics problems. The gravity can be considered as a constant force in most cases (see for dimensional analysis for the reasons).

Shear Stress and Surface Tension as Driving Force

Fig. 8.14 Kerosene lamp.

If the fluid was solid material, pulling the side will pull all the material. In fluid (mostly liquid) shear stress pulling side (surface) will have limited effect and yet sometime is significant and more rarely dominate. Consider, for example, the case shown in Figure 8.14. The shear stress carry the material as if part of the material was a solid material. For example, in the kerosene lamp the burning occurs at the surface of the lamp top and the liquid is at the bottom. The liquid does not move up due the gravity (actually it is against the gravity) but because the surface tension.

Fig. 8.15 Flow in a kendle with a surfece tension gradient.

The physical conditions in Figure 8.14 are used to idealize the flow around an inner rode to understand how to apply the surface tension to the boundary conditions. The fluid surrounds the rode and flows upwards. In that case, the velocity at the surface of the inner rode is zero. The velocity at the outer surface is unknown. The boundary condition at outer surface given by a jump of the shear stress. The outer diameter is depends on the surface tension (the larger surface tension the smaller the liquid diameter). The surface tension is a function of the temperature therefore the gradient in surface tension is result of temperature gradient. In this book, this effect is not discussed. However, somewhere downstream the temperature gradient is insignificant. Even in that case, the surface tension gradient remains. It can be noticed that, under the assumption presented here, there are two principal radii of the flow. One radius toward the center of the rode while the other radius is infinite (approximatly). In that case, the contribution due to the curvature is zero in the direction of the flow (see Figure 8.15). The only (almost) propelling source of the flow is the surface gradient ($\dfrac{\partial \sigma}{\partial n}$).

8.7 Examples for Differential Equation (Navier-Stokes)

Examples of an one-dimensional flow driven by the shear stress and pressure are presented. For further enhance the understanding some of the derivations are repeated. First, example dealing with one phase are present. Later, examples with two phase are presented.

Fig. 8.16 Flow between two plates, top plate is moving at speed of $U_\ll$ to the right (as positive). The control volume shown in darker colors.

Example 8.6

Incompressible liquid flows between two infinite plates from the left to the right (as shown in Figure 8.16). The distance between the plates is $\ell$. The static pressure per length is given as $\Delta P$ . The upper surface is moving in velocity, $U_{\ell}$ (The right side is defined as positive).

Solution

In this example, the mass conservation yields \begin{align} \label{dif:eq:2dmass} \overbrace{\dfrac{d}{dt} \int_{cv} \rho dV}^{=0} = - \int_{cv} \rho \,U_{rn} dA = 0 \end{align} The momentum is not accumulated (steady state and constant density). Further because no change of the momentum thus \begin{align} \label{diff:eq:whatDiffMass} \int_A \rho\,U_x\,U_{rn} dA = 0 \end{align} Thus, the flow in and the flow out are equal. It can concluded that the velocity in and out are the same (for constant density). The momentum conservation leads \begin{align} \label{dif:eq:2dmom} - \int_{cv} \pmb{P} dA + \int_{cv} \boldsymbol{\tau}_{xy} dA = 0 \end{align} The reaction of the shear stress on the lower surface of control volume based on Newtonian fluid is \begin{align} \label{dif:eq:2dlowerTau} \boldsymbol{\tau}_{xy} = - \mu \dfrac{dU}{dy} \end{align} On the upper surface is different by Taylor explanation as \begin{align} \label{dif:eq:eq:2dupperTau} \boldsymbol{\tau}_{xy} = \mu \left( \dfrac{dU}{dy} + \dfrac{d^2U}{dy^2}\,dy + \overbrace{\dfrac{d^3U}{dy^3}\,dy^2 +\cdots}^{\cong 0} \right) \end{align} The net effect of these two will be difference between them \begin{align} \label{dif:eq:muDif} \mu \left( \dfrac{dU}{dy\dfrac{}{\dfrac{}{}}} + \dfrac{d^2U}{dy^2}\,dy\right) - \mu \dfrac{dU}{dy} \cong \mu \, \dfrac{d^2U}{dy^2} \, dy \end{align} The assumptions is that there is no pressure difference in the $z$ direction. The only difference in the pressure is in the $x$ direction and thus \begin{align} \label{dif:eq:Pdif} P - \left( P + \dfrac{dP}{dx\dfrac{}{}} \,dx \right) = -\dfrac{dP}{dx} \,dx \end{align} A discussion why $\dfrac{\partial P}{\partial y} \sim 0$ will be presented later. The momentum equation in the $x$ direction (or from equation \eqref{dif:eq:momEqx}) results (without gravity effects) in \begin{align} \label{dif:eq:momX} -\dfrac{dP}{dx} = \mu \, \dfrac{d^2U}{dy^2} \end{align}

Fig. 8.17 One dimensional flow with a shear between two plates when $\Psi$ change value between -1.75 green line to 3 the blue line.

Equation \eqref{dif:eq:momX} was constructed under several assumptions which include the direction of the flow, Newtonian fluid. No assumption was imposed on the pressure distribution. Equation \eqref{dif:eq:momX} is a partial differential equation but can be treated as ordinary differential equation in the $z$ direction of the pressure difference is uniform. In that case, the left hand side is equal to constant. The standard'' boundary conditions is non–vanishing pressure gradient (that is the pressure exist) and velocity of the upper or lower surface or both. It is common to assume that the no slip'' condition on the boundaries condition . The boundaries conditions are \begin{align} \label{dif:eq:2Dbc} \displaystyle U_x(y=0) = 0 \\ \displaystyle U_x(y=\ell) = U_{\ell} \end{align} The solution of the ordinary'' differential equation \eqref{dif:eq:momX} after the integration becomes \begin{align} \label{dif:eq:2DGsolution} U_x = - \dfrac{1}{2} \dfrac{dP}{dx} \,y^2 + c_2\, y + c_3 \end{align} Applying the boundary conditions, equation \eqref{dif:eq:2Dbc}/ results in \begin{align} \label{dif:eq:2Dsolution} U_x (y) = \dfrac{y}{\ell} \left( \overbrace{\dfrac{\ell^2}{U_0\,2\mu} \dfrac{dP}{dx}}^{=\Psi} \left( 1 - \dfrac{y}{\ell} \right) \right) + \dfrac{y}{\ell} \end{align} For the case where the pressure gradient is zero the velocity is linear as was discussed earlier in Chapter 1 (see Figure \ref{dif:fig:oneDFlow}). However, if the plates or the boundary conditions do not move the solution is \begin{align} \label{dif:eq:2DsolutionNoPlate} U_x (y) = \left( \dfrac{\ell^2}{U_0\,2\mu\dfrac{}{}} \dfrac{dP}{dx} \left( 1 - \dfrac{y}{\ell} \right) \right) + \dfrac{y}{\ell} \end{align} What happen when $\dfrac{\partial P}{\partial y} \sim 0$?

Fig. 8.18 The control volume of liquid element in cylindrical coordinates.

Cylindrical Coordinates

Similarly the problem of one dimensional flow can be constructed for cylindrical coordinates. The problem is still one dimensional because the flow velocity is a function of (only) radius. This flow referred as Poiseuille flow after Jean Louis Poiseuille a French Physician who investigated blood flow in veins. Thus, Poiseuille studied the flow in a small diameters (he was not familiar with the concept of Reynolds numbers). Rederivation are carried out for a short cut. The momentum equation for the control volume depicted in the Figure 8.18a is \begin{align} \label{dif:eq:RdimMOM} - \int \pmb{P}\,dA + \int \boldsymbol{\tau} dA = \int \rho\,U_z\,U_{rn}\,dA \end{align} The shear stress in the front and back surfaces do no act in the $z$ direction. The shear stress on the circumferential part small dark blue shown in Figure 8.18a is \begin{align} \label{dif:eq:RdimtauR} \int \boldsymbol{\tau} \,dA = \mu \dfrac{dU_z}{dr} \overbrace{2\,\pi\,r\,dz}^{dA} \end{align} The pressure integral is \begin{align} \label{dif:eq:RdimP} \int \pmb{P}\,dA = \left( P_{z_dz} - P_{z} \right) \,\pi\,r^2 = \left( P_{z} + \dfrac{\partial P}{\partial z\dfrac{}{}} dz - P_z \right) \,\pi\,r^2 = \dfrac{\partial P}{\partial z} dz \,\pi\,r^2 \end{align} The last term is \begin{align*} \begin{array}{rl} \displaystyle \int \rho\,U_z\,U_{rn}\,dA = \rho\, \int U_z\,U_{rn}\,dA = \\ \displaystyle \rho \left( \dfrac{}{\dfrac{}{}} \int_{z+dz} \dfrac{}{} {U_{z+dz}}^2 dA \right. -\displaystyle \left. \dfrac{}{\dfrac{}{}} \int_{z} {U_{z}}^2 dA \right) = \displaystyle \rho \int_{z}\left( {U_{z+dz}}^2 - {U_{z}}^2 \right) dA \end{array} \end{align*} The term ${U_{z+dz}}^2 - {U_{z}}^2$ is zero because ${U_{z+dz}} = {U_{z}}$ because mass conservation conservation for any element. Hence, the last term is \begin{align} \label{dif:eq:RdimMOMlast} \int \rho\,U_z\,U_{rn}\,dA = 0 \end{align} Substituting equation \eqref{dif:eq:RdimtauR} and \eqref{dif:eq:RdimP} into equation \eqref{dif:eq:RdimMOM} results in \begin{align} \label{dif:eq:RdimMOMt0} \mu \dfrac{dU_z}{dr} 2\,\cancel{\pi}\,\cancel{r}\,\cancel{dz} = - \dfrac{\partial P}{\partial z} \cancel{dz} \,\cancel{\pi}\,r^{\cancel{2}} \end{align} Which shrinks to \begin{align} \label{dif:eq:RdimMOMt} \dfrac{2\, \mu}{r}\, \dfrac{dU_z}{dr} = - \dfrac{\partial P}{\partial z} \end{align} Equation \eqref{dif:eq:RdimMOMt} is a first order differential equation for which only one boundary condition is needed. The no slip'' condition is assumed \begin{align} \label{dif:eq:RdimNoSlip} U_z(r=R) = 0 \end{align} Where $R$ is the outer radius of pipe or cylinder. Integrating equation \eqref{dif:eq:RdimMOMt} results in \begin{align} \label{dif:eq:RdimNoSlipInt} U_z = - \dfrac{1}{\mu} \dfrac{\partial P}{\partial z}\,r^2 + c_1 \end{align} was eliminated due to the smart short cut. The integration constant obtained via the application of the boundary condition which is \begin{align} \label{dif:eq:RdimapplyBC1} c_1 = - \dfrac{1}{\mu} \dfrac{\partial P}{\partial z}\,R^2 \end{align} The solution is \begin{align} \label{dif:eq:RdimSsolution} U_z = \dfrac{1}{\mu} \dfrac{\partial P}{\partial z}\,R^2 \left(1 - \left( \dfrac{r}{R\dfrac{}{}} \right)^2 \right) \end{align} While the above analysis provides a solution, it has several deficiencies which include the ability to incorporate different boundary conditions such as flow between concentering cylinders.

Example 8.7

Fig. 8.19 Liquid flow between concentric cylinders for example .

A liquid with a constant density is flowing between concentering cylinders as shown in Figure . Assume that the velocity at the surface of the cylinders is zero calculate the velocity profile. Build the velocity profile when the flow is one directional and viscosity is Newtonian. Calculate the flow rate for a given pressure gradient.

Solution

After the previous example, the appropriate version of the Navier–Stokes equation will be used. The situation is best suitable to solved in cylindrical coordinates. One of the solution of this problems is one dimensional solution. In fact there is no physical reason why the flow should be only one dimensional. However, it is possible to satisfy the boundary conditions. It turn out that the simple'' solution is the first mode that appear in reality. In this solution will be discussing the flow first mode. For this mode, the flow is assumed to be one dimensional. That is, the velocity isn't a function of the angle, or $z$ coordinate. Thus only equation in $z$ coordinate is needed. It can be noticed that this case is steady state and also the acceleration (convective acceleration) is zero \begin{align} \label{CinC:eq:internalMom} \rho \left(\overbrace{\dfrac{\partial U_z}{\partial t}}^{\neq f(t)} + \overbrace{U_r}^{= 0} \dfrac{\partial U_z}{\partial r} + \overbrace{\dfrac{U_{\phi}}{r}}^{=0} \overbrace{\dfrac{\partial U_z}{\partial \phi}} ^{U_z\neq f(\phi)} + U_z \overbrace{\dfrac{\partial U_z}{\partial z}}^{=0}\right) = 0 \end{align} The steady state governing equation then becomes \begin{align} \label{dif:eq:govConcentricCylinders} \rho \left( \cancel{0} \right) = 0 = -\dfrac{\partial P}{\partial z} + \mu \left(\dfrac{1}{r}\dfrac{\partial}{\partial r}\left(r \dfrac{\partial U_z}{\partial r\dfrac{}{}}\right) + \overbrace{\cdots}^{=0}\right) + \cancel{\rho\, g_z} \end{align} The PDE above () required boundary conditions which are \begin{align*} U_z\,(r= r_i) = &0 \\ U_z\,(r= r_o) = &0 \end{align*} Integrating equation \eqref{dif:eq:govConcentricCylinders} once results in \begin{align} \label{CinC:eq:In1} r \,\dfrac{\partial U_z}{\partial r} = \dfrac{1}{2\,\mu} \dfrac{\partial P}{\partial z} r^2 + c_1 \end{align} Dividing equation \eqref{CinC:eq:In1} and integrating results for the second times results \begin{align} \label{CinC:eq:In1a} \dfrac{\partial U_z}{\partial r} = \dfrac{1}{2\,\mu} \dfrac{\partial P}{\partial z}\, r + \dfrac{c_1}{r} \end{align} Integration of equation \eqref{CinC:eq:In1a} results in \begin{align} \label{CinC:In2} U_z = \dfrac{1}{4\,\mu} \dfrac{\partial P}{\partial z} r^2 + c_1\, \ln {r} + c_2 \end{align} Applying the first boundary condition results in \begin{align} \label{CinC:bc1} 0 = \dfrac{1}{4\,\mu} \dfrac{\partial P}{\partial z} {r_i}^2 + c_1\, \ln {r_i} + c_2 \end{align} applying the second boundary condition yields \begin{align} \label{CinC:bc2} 0 = \dfrac{1}{4\,\mu} \dfrac{\partial P}{\partial z} {r_o}^2 + c_1\, \ln {r_o} + c_2 \end{align} The solution is \begin{align*} c_1 = \dfrac{1}{4\,\mu} \, \ln\left(\dfrac{r_o}{r_i}\right)\, \dfrac{\partial P}{dz \dfrac{}{}}\left( {r_o}^{2} - {r_i}^{2} \right) \\ c_2 = \dfrac{1}{4\,\mu} \, \ln\left(\dfrac{r_o}{r_i\dfrac{}{}}\right)\, \dfrac{\partial P}{dz }\left( \ln(r_i)\,{r_o}^{2} - \ln(r_o)\, {r_i}^{2} \right) \end{align*} The solution is when substituting the constants into equation \eqref{CinC:In2} results in \begin{align} \label{CinC:solution} \begin{array}{rcl} \displaystyle U_z(r) = \dfrac{1}{4\,\mu} \dfrac{\partial P}{\partial z} {r}^2 + \dfrac{1}{4\,\mu} \, \ln\left(\dfrac{r_o}{r_i}\right)\, \dfrac{\partial P}{dz }\left( {r_o}^{2} - {r_i}^{2} \right) \ln \, r \\ + \dfrac{1}{4\,\mu} \, \ln\left(\dfrac{r_o}{r_i}\right)\, \dfrac{\partial P}{dz }\left( \ln(r_i)\,{r_o}^{2} - \ln(r_o)\, {r_i}^{2} \right) \end{array} \end{align} The flow rate is then \begin{align} \label{CinC:FlowRateg} Q = \int_{r_i}^{r_o} U_z(r) dA \end{align} Or substituting equation \eqref{CinC:solution} into equation \eqref{CinC:FlowRateg} transformed into \begin{align} \label{CinC:Qsol} \begin{array}{rcl} \displaystyle Q = \int_A \left[ \dfrac{1}{4\,\mu} \dfrac{\partial P}{\partial z} {r}^2\right. + \dfrac{1}{4\,\mu} \, \ln\left(\dfrac{r_o}{r_i}\right)\, \dfrac{\partial P}{dz }\left( {r_o}^{2} - {r_i}^{2} \right) \ln \, r \\ + \dfrac{1}{4\,\mu} \, \ln\left(\dfrac{r_o}{r_i}\right)\, \left.\dfrac{\partial P}{dz }\left( \ln(r_i)\,{r_o}^{2} - \ln(r_o)\, {r_i}^{2} \right) \right] dA \end{array} \end{align} A finite integration of the last term in the integrand results in zero because it is constant. The integration of the rest is \begin{align} \label{CinC:Qsol1} \begin{array}{rcl} \displaystyle Q = \left[ \dfrac{1}{4\,\mu\dfrac{}{}} \dfrac{\partial P}{\partial z} \right] \int_{r_i}^{r_o} \left[ {r}^2 + \, \ln\left(\dfrac{r_o}{r_i}\right) \left( {r_o}^{2} - {r_i}^{2} \right) \ln \, r \right] 2\,\pi\,r\,dr \end{array} \end{align} The first integration of the first part of the second square bracket, ($r^3$), is $1/4\left({r_o}^4 - {r_i}^4 \right)$. The second part, of the second square bracket, $\left(-a\times r \ln \,r \right)$ can be done by parts to be \begin{align*} a\,\left(\dfrac{{r}^{2}}{4}-\dfrac{{r}^{2}\,\mathrm{log}\left( r\right) }{2\dfrac{}{}} \right) \end{align*} Applying all these techniques'' to equation \eqref{CinC:Qsol1} results in \begin{multline*} Q = \left[ \dfrac{\pi}{2\,\mu\dfrac{}{}} \dfrac{\partial P}{\partial z} \right] \left[ \left( \dfrac{{r_o}^{4}}{4 \dfrac{}{}}-\dfrac{{r_i}^{4}}{4} \right) \right. + \\ \displaystyle \ln\left(\dfrac{r_o}{r_i}\right) \left. \left( {r_o}^{2} - {r_i}^{2} \right) \left( \dfrac{{r_o}^{2}\,\mathrm{ln}\left( r_o\right)}{2} -\dfrac{{r_o}^{2}}{4}-\dfrac{{r_i}^{2}\,\mathrm{ln}\left( r_i\right) }{2\dfrac{}{}}+\dfrac{{r_i}^{2}}{4}\right) \right] \end{multline*} The averaged velocity is obtained by dividing flow rate by the area $Q/A$. \begin{align} \label{CinC:Uave} U_{ave} = \dfrac{Q}{ \pi\left( {r_o}^{2} - {r_i}^{2} \right) } \end{align} in which the identity of $(a^4-b^4)/(a^2-b^2)$ is $b^2+a^2$ and hence \begin{multline*} \displaystyle U_{ave} = \left[ \dfrac{1}{2\,\mu} \dfrac{\partial P}{\partial z\dfrac{}{}} \right] \left[ \left( \dfrac{{r_o}^{2}}{4\dfrac{}{}}+\dfrac{{r_i}^{2}}{4} \right) \right. + \\ \displaystyle \ln\left(\dfrac{r_o}{r_i\dfrac{}{}}\right) \left. \left( \dfrac{{r_o}^{2}\,\mathrm{ln}\left( r_o\right)}{2} -\dfrac{{r_o}^{2}}{4}-\dfrac{{r_i}^{2}\,\mathrm{ln}\left( r_i\right) }{2\dfrac{}{}}+\dfrac{{r_i}^{2}}{4}\right) \right] \end{multline*}

Example 8.8

For the conteraic veloicty profile, at what radius the maximum velocity obtained. Draw the maximum velocity location as a funciton of the ratio $r_{i}/r_o$.

The next example deals with the gravity as body force in two dimensional flow. This problem study by Nusselt which developed the basics equations. This problem is related to many industrial process and is fundamental in understanding many industrial processes. Furthermore, this analysis is a building bloc for heat and mass transfer understanding .

Example 8.9

Fig. 8.20 Mass flow due to temperature difference for example

In many situations in nature and many industrial processes liquid flows downstream on inclined plate at $\theta$ as shown in Figure 8.20. For this example, assume that the gas density is zero (located outside the liquid domain). Assume that scale'' is large enough so that the no slip'' condition prevail at the plate (bottom). For simplicity, assume that the flow is two dimensional.} Assume that the flow obtains a steady state after some length (and the acceleration vanished). The dominate force is the gravity. Write the governing equations for this situation. Calculate the velocity profile. Assume that the flow is one dimensional in the $x$ direction.

Solution

This problem is suitable to Cartesian coordinates in which $x$ coordinate is pointed in the flow direction and $y$ perpendicular to flow direction (depicted in Figure 8.20). For this system, the gravity in the $x$ direction is $g\sin\theta$ while the direction of $y$ the gravity is $g\cos\theta$. The governing in the $x$ direction is \begin{align} \label{gravity:gov} \begin{array} {rll} \rho \left(\overbrace{\dfrac{\partial U_x}{\partial t}}^{\neq f(t)} + \right. & U_x\overbrace{\dfrac{\partial U_x}{\partial x}}^{=0} + \overbrace{U_y}^{=0} \dfrac{\partial U_x}{\partial y} + \left. \overbrace{U_z}^{-0} \dfrac{\partial U_x}{\partial z}\right) = &\\ &-\overbrace{\dfrac{\partial P}{\partial x}}^{\sim 0} + \mu \left(\overbrace{\dfrac{\partial^2 U_x}{\partial x^2}}^{=0} + \dfrac{\partial^2 U_x}{\partial y^2} + \overbrace{\dfrac{\partial^2 U_x}{\partial z^2}}^{=0}\right) + \rho \overbrace{g_x}^{g\sin\theta} \end{array} \end{align} The first term of the acceleration is zero because the flow is in a steady state. The first term of the convective acceleration is zero under the assumption of this example flow is fully developed and hence not a function of $x$ (nothing to be improved''). The second and the third terms in the convective acceleration are zero because the velocity at that direction is zero ($U_y=U_z=0$). The pressure is almost constant along the $x$ coordinate. As it will be shown later, the pressure loss in the gas phase (mostly air) is negligible. Hence the pressure at the gas phase is almost constant hence the pressure at the interface in the liquid is constant. The surface has no curvature and hence the pressure at liquid side similar to the gas phase and the only change in liquid is in the $y$ direction. Fully developed flow means that the first term of the velocity Laplacian is zero ($\dfrac{\partial U_x}{\partial x} quiv 0$). The last term of the velocity Laplacian is zero because no velocity in the $z$ direction. Thus, equation \eqref{gravity:gov} is reduced to \begin{align} \label{gravity:govRed} 0 = \mu \, \dfrac{\partial^2 U_x}{\partial y^2} + \rho \,g\sin\theta \end{align} With boundary condition of no slip'' at the bottom because the large scale and steady state \begin{align} \label{gravity:bcBotton} U_x (y=0) =0 \end{align} The boundary at the interface is simplified to be \begin{align} \label{gravity:bcSurface1} \left.\dfrac{\partial U_x}{\partial y} \right|_{y=0} = \tau_{air} \left( \sim 0 \right) \end{align} If there is additional requirement, such a specific velocity at the surface, the governing equation can not be sufficient from the mathematical point of view. Integration of equation \eqref{gravity:govRed} yields \begin{align} \label{gravity:govInt} \dfrac{\partial U_x}{\partial y} = \dfrac{\rho}{\mu} \,g\sin\theta\,y + c_1 \end{align} The integration constant can be obtain by applying the condition \eqref{gravity:bcSurface1} \begin{align} \label{gravity:govApl1} \tau_{air} = \mu \left. \dfrac{\partial U_x}{\partial y}\right|_h = - \rho \,g\sin\theta\,\overbrace{h}^{y} + c_1 \, \mu \end{align} Solving for $c_1$ results in \begin{align} \label{gravity:c1} c_1 = \dfrac{\tau_{air} }{ \mu} + \dfrac{1}{\underbrace{\nu}_{\dfrac{\mu}{\rho}}} \, g\,\sin\theta\,h \end{align} The second integration applying the second boundary condition yields $c_2=0$ results in \begin{align} \label{gravity:govInt2} U_x = \dfrac{g\,\sin\theta}{\nu} \left( 2\,y\,h- y^2 \right) - \dfrac{\tau_{air}}{\mu} \end{align} When the shear stress caused by the air is neglected, the velocity profile is \begin{align} \label{gravity:solNoShear} U_x = \dfrac{g\,\sin\theta}{\nu} \left( 2\,h\,y- y^2 \right) \end{align} The flow rate per unit width is \begin{align} \label{gravity:flowRateIn} \dfrac{Q}{W} = \int_{A} U_x dA = \int_0^h \left( \dfrac{g\,\sin\theta}{\nu} \left( 2\,h\,y - y^2 \right) - \dfrac{\tau_{air}}{\mu}\, \right) dy \end{align} Where $W$ here is the width into the page of the flow. Which results in \begin{align} \label{gravity:flowRate1} \dfrac{Q}{W} = \dfrac{g\,\sin\theta}{\nu} \dfrac{2\,h^3}{3} - \dfrac{\tau_{air}\,h}{\mu} \end{align} The average velocity is then \begin{align} \label{gravity:flowRate} \overline{U_x} = \dfrac{\dfrac{Q}{W}}{h} = \dfrac{g\,\sin\theta}{\nu} \dfrac{2\,h^2}{3} - \dfrac{\tau_{air}}{\mu} \end{align} Note the shear stress at the interface can be positive or negative and hence can increase or decrease the flow rate and the averaged velocity.

In the following following example the issue of driving force of the flow through curved interface is examined. The flow in the kerosene lamp is depends on the surface tension. The flow surface is curved and thus pressure is not equal on both sides of the interface.

Example 8.10

A simplified flow version the kerosene lump is of liquid moving up on a solid core. Assume that radius of the liquid and solid core are given and the flow is at steady state. Calculate the minimum shear stress that required to operate the lump (alternatively, the maximum height).

8.7.1 Interfacial Instability

Fig. 8.? Flow of liquid in partially filled duct.

In Example 8.9 no requirement was made as for the velocity at the interface (the upper boundary). The vanishing shear stress at the interface was the only requirement was applied. If the air is considered two governing equations must be solved one for the air (gas) phase and one for water (liquid) phase. Two boundary conditions must be satisfied at the interface. For the liquid, the boundary condition of no slip'' at the bottom surface of liquid must be satisfied. Thus, there is total of three boundary conditions to be satisfied. The solution to the differential governing equations provides only two constants. The second domain (the gas phase) provides another equation with two constants but again three boundary conditions need to satisfied. However, two of the boundary conditions for these equations are the identical and thus the six boundary conditions are really only 4 boundary conditions. The governing equation solution for the gas phase ($h \geq y \geq a\,h$) is \begin{align} \label{bmI:eq:gas} {U_x}_g = \dfrac{g\sin\theta}{2\,\nu_g} \,y^2 + c_1\,y + c_2 \end{align} Note, the constants $c_1$ and $c_2$ are dimensional which mean that they have physical units ($c_1 \longrightarrow \left[1/sec\right]$ The governing equation in the liquid phase ($0 \geq y \geq h$) is \begin{align} \label{bmI:eq:liquid} {U_x}_{\ell} = \dfrac{g\sin\theta}{2\,\nu_{\ell}} \,y^2 + c_3\,y + c_4 \end{align} The gas velocity at the upper interface is vanished thus \begin{align} \label{bmI:eq:gbc1} {U_x}_g\left[ (1+a)\,h\right] = 0 \end{align} At the interface the no slip'' condition is regularly applied and thus \begin{align} \label{bmI:eq:gbc2} {U_x}_g(h) = {U_x}_{\ell} (h) \end{align} Also at the interface (a straight surface), the shear stress must be continuous \begin{align} \label{bmI:eq:gbc3} \mu_{g} \dfrac{\partial {U_x}_g}{\partial y} = \mu_{\ell} \dfrac{\partial {U_x}_{\ell}}{\partial y} \end{align} Assuming no slip'' for the liquid at the bottom boundary as \begin{align} \label{bmI:eq:lbc1} {U_x}_{\ell} (0) = 0 \end{align} The boundary condition \eqref{bmI:eq:gbc1} results in \begin{align} \label{bmI:eq:gbcs1} 0 = \dfrac{g\,\sin\theta}{2\,\nu_g}\,h^2\,(1+a)^2 + c_1 \,h\,(1+a) + c_2 \end{align} The same can be said for boundary condition \eqref{bmI:eq:lbc1} which leads \begin{align} \label{bmI:eq:gbcs1c} c_4 = 0 \end{align} Applying equation \eqref{bmI:eq:gbc3} yields \begin{align} \label{bmI:eq:gbcs3} \overbrace{\dfrac{\mu_g}{\nu_g}}^{\rho_g} {g\,\sin\theta}\, h + c_1\,\mu_g= \overbrace{\dfrac{\mu_{\ell}}{\nu_{\ell}}}^{\rho_{\ell}} {g\,\sin\theta}\, h + c_3\,\mu_{\ell} \end{align} Combining boundary conditions equationqref{bmI:eq:gbc2} with \eqref{bmI:eq:gbcs1} results in \begin{align} \label{bmI:eq:gbcs2} \dfrac{g\sin\theta}{2\,\nu_{g}} \,h^2 + c_1\,h + c_2 = \dfrac{g\sin\theta}{2\,\nu_{\ell}} \,h^2 + c_3\,h \end{align}

The solution of equation \eqref{bmI:eq:gbcs1}, \eqref{bmI:eq:gbcs3} and \eqref{bmI:eq:gbcs2} is obtained by computer algebra (see in the code) to be \begin{multline} \label{bmI:eq:gSol} \displaystyle c_1= -\dfrac{\sin\theta\,\left( g\,h\,\rho_{g}\,\left( 2\,\rho_{g}\,\nu_{\ell}\,\rho_{\ell}+1\right) +a\,g\,h\,\nu_{\ell}\right) } {\rho_{g}\,\left( 2\,a\,\nu_{\ell}+2\,\nu_{\ell}\right) } \\ \displaystyle c_2= \dfrac{\sin\theta\,\left( g\,{h}^{2}\,\rho_{g}\,\left( 2\,\rho_{g}\,\nu_{\ell}\,\rho_{\ell}+1\right) -g\,{h}^{2}\,\nu_{\ell}\right) }{2\,\rho_{g}\,\nu_{\ell}} \\ \displaystyle c_3= \dfrac{\sin\theta\,\left( g\,h\,\rho_{g}\,\left( 2\,a\,\rho_{g}\,\nu_{\ell}\,\rho_{\ell}-1\right) -a\,g\,h\,\nu_{\ell}\right) }{\rho_{g}\,\left( 2\,a\,\nu_{\ell}+2\,\nu_{\ell}\right) } \end{multline}

When solving this kinds of mathematical problem the engineers reduce it to minimum amount of parameters to reduce the labor involve. So equation \eqref{bmI:eq:gbcs1} transformed by some simple rearrangement to be \begin{align} \label{bmI:eq:gbcsd1} \left(1+a\right)^{2} = \overbrace{\dfrac{2\,\nu_g\,c_1}{g\,h\,\sin\theta}}^{C_1} + \overbrace{\dfrac{2\,c_2 \,\nu_g}{g\,h^2\,\sin\theta}}^{C_2} \end{align} And equation \eqref{bmI:eq:gbcs3} \begin{align} \label{bmI:eq:gbcsd3} 1 + \overbrace{\dfrac{\nu_g\, c_1}{g\,h\,\sin\theta}}^{\dfrac{1}{2}\,C_1} = \dfrac{\rho_{\ell}} {\rho_g} + \overbrace{\dfrac{\mu_{\ell}\,\nu_g\, c_3}{\mu_{g}\,g\,h\,\sin\theta}} ^{\dfrac{1}{2}\dfrac{\mu_{\ell}}{\mu_g}\,C_3 } \end{align} and equation \eqref{bmI:eq:gbcs2} \begin{align} \label{bmI:eq:gbcsd2} 1 +\dfrac{2\,\nu_g\,\cancel{h}\,c_1}{h^{\cancel{2}}\,g\,\sin\theta} + \dfrac{2\,\nu_{g}\,c_2} {h^2\,g\sin\theta} = \dfrac{\nu_g}{\nu_{\ell}} + \dfrac{2\,\nu_{g}\,\cancel{h}\,c_3} {g\,h^{\cancel{2}}\,\sin\theta} \end{align} Or rearranging equation \eqref{bmI:eq:gbcsd2} \begin{align} \label{bmI:eq:gbcsd2a} \dfrac{\nu_g}{\nu_{\ell}} - 1 = \overbrace{\dfrac{2\,\nu_g\,c_1}{h\,g\,\sin\theta}}^{C_1} + \overbrace{\dfrac{2\,\nu_{g}\,c_2} {h^2\,g\sin\theta}}^{C_2} - \overbrace{\dfrac{2\,\nu_{g}\,\,c_3} {g\,h\,\sin\theta} }^{C_3} \end{align} This presentation provide similarity and it will be shown in the Dimensional analysis chapter better physical understanding of the situation. Equation \eqref{bmI:eq:gbcsd1} can be written as \begin{align} \label{bmI:eq:gbcsdd1} \left(1+a\right)^2 = C_1 + C_2 \end{align} Further rearranging equation \eqref{bmI:eq:gbcsd3} \begin{align} \label{bmI:eq:gbcsdf3} \dfrac{\rho_{\ell}} {\rho_g} -1 = \dfrac{ C_1}{2} - \dfrac{\mu_{\ell}}{\mu_g}\, \dfrac{C_3}{2} \end{align} and equation \eqref{bmI:eq:gbcsd2a} \begin{align} \label{bmI:eq:gbcsdd2} \dfrac{\nu_g}{\nu_{\ell}} - 1 = C_1 + C_2 -C_3 \end{align} This process that was shown here is referred as non–﻿dimensionalization . The ratio of the dynamics viscosity can be eliminated from equation \eqref{bmI:eq:gbcsdd2} to be \begin{align} \label{bmI:eq:gbcsdd2m} \dfrac{\mu_g}{\mu_{\ell}} \,\dfrac{\rho_{\ell}} {\rho_g} - 1 = C_1 + C_2 -C_3 \end{align} The set of equation can be solved for the any ratio of the density and dynamic viscosity. The solution for the constant is \begin{align} \label{bmI:eq:solC1} C_1 = \dfrac{\rho_{g}}{\rho_{\ell}}-2 +{a}^{2}+2\,a\dfrac{\mu_{g}}{\mu_{\ell}}+2\dfrac{\mu_{g}}{\mu_{\ell}} \end{align} \begin{align} \label{bmI:eq:solC2} C_2=\dfrac{-\dfrac{\mu_{g}}{\mu_{\ell}}\,\dfrac{\rho_{\ell}}{\rho_g}+ a\,\left( 2\,\dfrac{\mu_{g}}{\mu_{\ell}}-2\right) + 3\,\dfrac{\mu_{g}}{\mu_{\ell}}+{a}^{2}\,\left( \dfrac{\mu_{g}}{\mu_{\ell}}-1\right) -2}{\dfrac{\mu_{g}}{\mu_{\ell}}} \end{align} \begin{align} \label{bmI:eq:solC3} C_3=-\dfrac{\mu_{g}}{\mu_{\ell}}\,\dfrac{\rho_{\ell}}{\rho_g}+{a}^{2}+2\,a+2 \end{align} The two different fluids have flow have a solution as long as the distance is finite reasonable similar. What happen when the lighter fluid, mostly the gas, is infinite long. This is one of the source of the instability at the interface. The boundary conditions of flow with infinite depth is that flow at the interface is zero, flow at infinite is zero. The requirement of the shear stress in the infinite is zero as well. There is no way obtain one dimensional solution for such case and there is a component in the $y$ direction. Combining infinite size domain of one fluid with finite size on the other one side results in unstable interface.

Spherical Coordinates

Next: Dimensional1    Previous: Differential2