Category: Advanced Chemistry

Posted on

Change of basis

A change of basis is a method to convert vector coordinates with respect to one basis to coordinates with respect to another basis.

Consider a vector u with reference to orthogonal unit basis vectors x1, x2 and x3.

\textbf{\textit{u}}=px_1+qx_2+rx_3\; \; \; \; \; \; \; (1)

where p, q and r are the scalar components of the unit basis vectors x1, x2 and x3 respectively.

Clearly,

p=\textbf{\textit{u}}\cdot x_1\; \; \; \; q=\textbf{\textit{u}}\cdot x_2\; \; \; \; r=\textbf{\textit{u}}\cdot x_3\; \; \; \; (2)

If we define u with respect to another orthogonal reference frame with unit basis vectors x’1, x’2 and x’3 (see above diagram), we have:

\textbf{\textit{u}}=p' x'_1+q'x'_2+r'x'_3\; \; \; \; (3)

and

p'=\textbf{\textit{u}}\cdot x'_1\; \; \; \; q'=\textbf{\textit{u}}\cdot x'_2\; \; \; \; r'=\textbf{\textit{u}}\cdot x'_3\; \; \; \; (4)

To find the relationship between the two reference frames, we substitute eq1 in eq4 to give

p'=\left ( px_1+qx_2+rx_3 \right )\cdot x'_1=px_1\cdot x'_1+qx_2\cdot x'_1+rx_3\cdot x'_1\; \; \; \; \; \; \; (5)

q'=\left ( px_1+qx_2+rx_3 \right )\cdot x'_2=px_1\cdot x'_2+qx_2\cdot x'_2+rx_3\cdot x'_2\; \; \; \; \; \; \; (6)

r'=\left ( px_1+qx_2+rx_3 \right )\cdot x'_3=px_1\cdot x'_3+qx_2\cdot x'_3+rx_3\cdot x'_3\; \; \; \; \; \; \; (7)

Since the dot product of two vectors gives a scalar, we can write eq5, eq6 and eq7 as

p'=a_{11}p+a_{12}q+a_{13}r

q'=a_{21}p+a_{22}q+a_{23}r

r'=a_{31}p+a_{32}q+a_{33}r

which can be expressed in the following matrix equation:

\begin{pmatrix} p'\\ q'\\ r' \end{pmatrix}=\begin{pmatrix} a_{11} &a_{12} &a_{13} \\ a_{21} & a_{22} &a_{23} \\ a_{31} &a_{32} & a_{33} \end{pmatrix}\begin{pmatrix} p\\ q\\ r \end{pmatrix}\; \; \; \; \; \; \; (8)

where a_{ij}=x'_i\cdot x_j=\left | x'_i \right |\left | x_j \right |cos\theta_{x'_i,x_j}=cos\theta _{x'_i,x_j}.

The matrix containing a­ij is known as the change of basis matrix . We can rewrite eq8 as:

u'_i=\sum_{j=1,2,3}a_{ij}u_j\; \; \; \; \; \; \; (9)

where u'_1=p',u'_2=q',u'_3=r' and u_1=p,u_2=q,u_3=r.

Eq9 is often written as a short form by omitting the summation symbol:

u'_i=a_{ij}u_j\; \; \; \; \; \; \; (10)

For example, the change of basis from the orthogonal basis xj to the orthogonal basis  x’i as a rotation about the x3-axis, where x3 coincides with x’3 (x’3 = x3 and is perpendicular to the plane of the page), is depicted as:

To determine the change of basis matrix, we have

x'_1\cdot x_3=x'_3\cdot x_1=x'_2\cdot x_3=x'_3\cdot x_2=cos90^o=0

x'_3\cdot x_3=cos0^o=1

x'_1\cdot x_1=x'_2\cdot x_2=cos\theta

x'_1\cdot x_2=cos\left ( 90^o-\theta \right )=sin\theta

x'_2\cdot x_1=cos\left ( 90^o+\theta \right )=-sin\theta

Therefore,

a_{ij}=\begin{pmatrix} a_{11} & a_{12}&a_{13}\\ a_{21}& a_{22}&a_{23}\\ a_{31}& a_{32}& a_{33}\end{pmatrix}=\begin{pmatrix} x'_1\cdot x_1 & x'_1\cdot x_2&x'_1\cdot x_3\\ x'_2\cdot x_1&x'_2\cdot x_2 &x'_2\cdot x_3\\ x'_3\cdot x_1& x'_3\cdot x_2& x'_3\cdot x_3\end{pmatrix}=\begin{pmatrix} cos\theta &sin\theta &0 \\ -sin\theta & cos\theta &0 \\ 0 &0 &1 \end{pmatrix}\; \; \; \; \; \; \; (12)

Consider again the vector u with reference to a coordinate system that is defined by the orthogonal unit basis vectors x1, x2 and x3  (see diagram below).

The rotation of u with respect to this fixed coordinate system is expressed as:

Compared to the change of basis of u that was described earlier, the rotation matrix R is the inverse of the change of basis matrix (similarly, a reflection matrix is the inverse of the change of basis matrix by a reflection). In other words, the rotation of u by θ can be perceived a change of basis, with the coordinate system rotated by –θ. This implies that we can analyse a problem in two ways.

 

Next article: Tensor and tensor transformation
Content page of tensors
Content page of advanced chemistry
Main content page
Posted on

Constitutive relation: hypothesis 3

The third constitutive relation hypothesis states that the shear stress tensor is zero if the flow involves no shearing of the fluid-body.

A square matrix A is symmetric if aij = aji for all i and j, i.e. A = AT, e.g.

\begin{pmatrix} 1 &5 &-3 \\ 5 &8 &0 \\ -3 &0 &2 \end{pmatrix}

A square matrix is anti-symmetric (or skew symmetric) if aij = –aji for all i and j, i.e. –A = AT, e.g.

\begin{pmatrix} 0 &-7 &3 \\ 7 &0 &-4 \\ -3 &4 &0 \end{pmatrix}

If aij = –aji , the diagonal components of the matrix must all be zero.

In general, we can write a square matrix in the following form:

A=\frac{1}{2}\left ( A+A^T \right )+\frac{1}{2}\left ( A-A^T \right )=B+C

where B = (A + AT)/2, C = (A – AT)/2 and AT is the transpose of A.

Using the transpose identity of (X + Y)T = XT+ YT (which can be verified by substituting any matrix with any components into the equality), we can show that B is symmetric (B = BT) and C is anti-symmetric (C = –CT).

This means that a square matrix can be decomposed into a symmetric and an anti-symmetric component. From the previous section, \frac{\partial u_k}{\partial x_l} is a second-order tensor that can be represented by a 3×3 matrix. Therefore,

\frac{\partial u_k}{\partial x_l}=\left [ \frac{1}{2}\left ( \frac{\partial u_k}{\partial x_l}+\frac{\partial u_l}{\partial x_k}\right )+\frac{1}{2}\left ( \frac{\partial u_k}{\partial x_l}-\frac{\partial u_l}{\partial x_k}\right )\right ]\; \; \; \; \; \; \; (9)

where \frac{1}{2}\left ( \frac{\partial u_k}{\partial x_l}+\frac{\partial u_l}{\partial x_k}\right ) is symmetric and \frac{1}{2}\left ( \frac{\partial u_k}{\partial x_l}-\frac{\partial u_l}{\partial x_k}\right ) is anti-symmetric.

Substitute eq9 in eq8,

\tau_{ij}=\alpha_{ijkl}\left [ \frac{1}{2}\left ( \frac{\partial u_k}{\partial x_l}+\frac{\partial u_l}{\partial x_k}\right )+\frac{1}{2}\left ( \frac{\partial u_k}{\partial x_l}-\frac{\partial u_l}{\partial x_k}\right )\right ]\; \; \; \; \; \; \; (10)

Consider a two-dimensional body of fluid ABCD as depicted in the above diagram. The term \left ( \frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1}\right ) tends to shear the fluid. If we consider the rotation of the fluid in the clockwise direction, the term \left ( \frac{\partial u_1}{\partial x_2}-\frac{\partial u_2}{\partial x_1}\right ) tends to rotate the fluid. This means that \frac{\partial u_k}{\partial x_l} in eq9 can be non-zero even if the shear component \left ( \frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1}\right ) is zero, as the rotation component \left ( \frac{\partial u_1}{\partial x_2}-\frac{\partial u_2}{\partial x_1}\right ) is non-zero. However, hypothesis 3 states that the shear stress tensor τij in eq10 is zero if the flow involves no shearing of the fluid-body. To satisfy hypothesis 3, eq10 becomes:

\tau_{ij}=\beta_{ijkl}\: \frac{1}{2}\left ( \frac{\partial u_k}{\partial x_l}+\frac{\partial u_l}{\partial x_k}\right )\; \; \; \; \; \; \; (11)

βijkl , which is composed of eighty-one components, characterizes the properties of the fluid. Since  \frac{1}{2}\left ( \frac{\partial u_k}{\partial x_l}+\frac{\partial u_l}{\partial x_k}\right ) is symmetric, and βijkl  consists of scalar components, τij is also symmetric.

 

next article: hypothesis 4
Previous article: hypothesis 2
Content page of constitutive relation
Content page of advanced chemistry
Main content page
Posted on

Constitutive relation: overview

A constitutive relation is an equation expressing the relation between two physical quantities (usually one kinetic and the other kinematic) with reference to a material. For example, Hooke’s law F = -kx states that the force needed to compress a particular spring by a certain distance varies linearly with that distance.

For an incompressible, viscous Newtonian fluid, George Stokes proposed that the constitutive relation must satisfy the following hypotheses:

  1. The stress tensor of a fluid at rest is the pressure exerted by the fluid (hydrostatic).
  2. The shear stress tensor is a linear function of the deformation tensor (velocity gradient).
  3. The shear stress tensor is zero if the flow involves no shearing of the fluid-body.
  4. The fluid is isotropic.

 

next article: hypothesis 1
Content page of constitutive relation
Content page of advanced chemistry
Main content page
Posted on

Constitutive relation: hypothesis 1

The first hypothesis, formulating the constitutive relation for an incompressible and viscous Newtonian fluid, states that the stress tensor of a fluid at rest is the pressure exerted by the fluid .

Stress is the measure of the amount of forces exerted by neighbouring particles of a material (e.g. a body of fluid) on one another. In fluid dynamics, stress consists of a hydrostatic portion and a hydrodynamic portion:

Stress=hydrostatic+hydrodynamic\; \; \; \; \; \; \; (1)

The components of stress that are perpendicular to the material’s cross section are simply the pressure exerted by the fluid (hydrostatic), while those that are coplanar with the material’s cross section are called shear stress (hydrodynamic).

To quantify a component of stress, we need to indicate a magnitude and two directions, one for the surface (surface direction is given by the normal vector to the surface) and the other for the direction of force acting on that surface. We can combine the two directions using a 2-letter subscript, with the first letter indicating the surface direction and the second letter denoting the force direction. For example, τxy is the component of stress with a y-direction force acting on the x-direction surface (i.e., surface BDFH). To further distinguish surface BDFH from surface ACEG, we write the y-direction force on BDFH as τxy(x+dx) and that on ACEG as τxy(x). 

Therefore, stress can be mathematically represented by a matrix that has nine components. We call such a matrix, a second-order tensor σij . To differentiate perpendicular components from coplanar components, we denote the former by σij and the latter by τij. The stress tensor can be expressed as:

\begin{pmatrix} \sigma_{xx} &\tau_{xy} &\tau_{xz} \\ \tau_{yx} &\sigma_{yy} &\tau_{yz} \\ \tau_{zx} &\tau_{zy} &\sigma_{zz} \end{pmatrix}\; \; \; \; or\; \; \; \; \sigma_{ij}=\begin{pmatrix} \sigma_{11} &\sigma_{12} &\sigma_{13} \\ \sigma_{21} &\sigma_{22} &\sigma_{23} \\ \sigma_{31} &\sigma_{32} &\sigma_{33} \end{pmatrix}

 

Question

How do we indicate stress components on the remaining 3 surfaces, i.e. ACEG, CDFE and ABDC?

Answer

Let’s illustrate this with an example. As mentioned above, to distinguish the y-direction force on BDFH from the y-direction force on ACEG, we denote each component as a function of distance, writing the former as τxy(x+dx) and the latter as τxy(x).

 

Hypothesis 1 states that, when the fluid is at rest:

\sigma_{11}=\sigma_{22}=\sigma_{33}=-p\; \; \; \; \; \; \; (2)

By convention, pressure is given a negative sign while stress, which is tensile in nature, is afforded a positive sign. Combining eq1 and eq2,

\sigma_{ij}=-p\delta_{ij}+\tau_{ij}\; \; \; \; \; \; \; (3)

where τij = 0 when the fluid is at rest and δij is the Kronecker delta.

Eq3 in matrix form is,

\begin{pmatrix} \sigma_{11} &\sigma_{12} &\sigma_{13} \\ \sigma_{21} &\sigma_{22} &\sigma_{23} \\ \sigma_{31} &\sigma_{32} &\sigma_{33} \end{pmatrix} =-p\begin{pmatrix} 1&0 &0 \\ 0 &1 &0 \\ 0 &0 &1 \end{pmatrix} +\begin{pmatrix} \tau_{11} &\tau_{12} &\tau_{13} \\ \tau_{21} &\tau_{22} &\tau_{23} \\ \tau_{31} &\tau_{32} &\tau_{33} \end{pmatrix} \; \; \; \; \; \; \; (4)

 

next article: hypothesis 2
Previous article: Overview
Content page of constitutive relation
Content page of advanced chemistry
Main content page
Posted on

Derivation of Stokes’ law

The derivation of Stokes’ law involves a few steps. Firstly, we shall derive the expression  \frac{\partial p}{\partial r}=\frac{\mu}{r^2sin\theta}\frac{\partial }{\partial \theta}E^2\psi. Substituting eq33b in eq33, we have:

\nabla p=-\mu\nabla\times\left ( \nabla\times\textbf{\textit{u}}\right )\; \; \; \; \; \; \; (47)

Substituting eq35b in eq47 and working out the algebra, we get:

\nabla p=\mu\left ( \frac{1}{r^2sin\theta} \frac{\partial }{\partial \theta}E^2\psi\right )\hat{r}-\mu\left ( \frac{1}{rsin\theta}\frac{\partial }{\partial r}E^2\psi \right )\hat{\theta}

The gradient of a function p in spherical coordinates is \nabla p= \frac{\partial p}{\partial r}\hat{r}+\frac{1}{r}\frac{\partial p}{\partial \theta}\hat{\theta}+\frac{1}{rsin\theta}\frac{\partial p}{\partial \phi}\hat{\phi}. Since fluid flow is along the z-axis, \nabla p= \frac{\partial p}{\partial r}\hat{r}+\frac{1}{r}\frac{\partial p}{\partial \theta}\hat{\theta}, which when compared with the above equation:

\frac{\partial p}{\partial r}=\frac{\mu}{r^2sin\theta}\frac{\partial }{\partial \theta}E^2\psi\; \; \; \; and\; \; \; \; \frac{1}{r}\frac{\partial p}{\partial \theta}=-\frac{\mu}{rsin\theta}\frac{\partial }{\partial r}E^2\psi\; \; \; \;\; \; \; (48)

Next, we shall show that E^2\psi=\frac{3}{2}\textbf{\textit{u}}\frac{a}{r}sin^2\theta. From a previous article, E^2=\frac{\partial ^2}{\partial r^2}+\frac{sin\theta}{r^2}\frac{\partial }{\partial \theta}\left ( \frac{1}{sin\theta}\frac{\partial }{\partial \theta} \right ), and so,

E^2\psi=\frac{\partial ^2}{\partial r^2}\psi+\frac{sin\theta}{r^2}\frac{\partial }{\partial \theta}\left ( \frac{1}{sin\theta}\frac{\partial }{\partial \theta} \right )\psi\; \; \; \; \; \; \; (49)

Substituting eq44 in eq49 and working out the algebra, we have:

E^2\psi=\frac{3}{2}\textbf{\textit{u}}\frac{a}{r}sin^2\theta\; \; \; \; \; \; \; (50)

Substituting eq50 in the first equation of eq48 and integrating, we have:

\int_{p}^{p_\infty }dp=\int_{r}^{r_\infty }\left [ \frac{\mu}{r^2sin\theta}\frac{\partial }{\partial \theta}\left ( \frac{3}{2} \textbf{\textit{u}}\frac{a}{r}sin^2\theta\right ) \right ]dr

p=p_\infty -\frac{3\mu\textbf{\textit{u}}a}{2r^2}cos\theta\; \; \; \; \; \; \; (51)

With reference to the diagram above, the stress τ on the sphere in the z-direction is given by the sum of the projections of τrr and τrθ on the z-axis:

\tau=\tau_{rr}cos\theta+\left [ -\tau_{r\theta}cos\left ( 90^o-\theta \right ) \right ]

\tau=\tau_{rr}cos\theta-\tau_{r\theta}sin\theta\; \; \; \; \; \; \; (52)

From the articles on constitutive relation, \sigma_{ij}=-p\delta_{ij}+\mu\left ( \frac{\partial u_i}{\partial x_j} +\frac{\partial u_j}{\partial x_i}\right ). So,

\tau_{rr}=-p+2\mu\frac{\partial u_r}{\partial r}\; \; \; \; \; \; \; (53)

\tau_{r\theta}=\mu\left ( r\frac{\partial }{\partial r}\frac{u_\theta}{r}+\frac{1}{r}\frac{\partial u_r}{\partial \theta} \right )\; \; \; \; \; \; \; (53a)

 

Question

Show that \tau_{r\theta}=\mu\left ( r\frac{\partial }{\partial r}\frac{u_\theta}{r}+\frac{1}{r}\frac{\partial u_r}{\partial \theta} \right ).

Answer

We need to convert the term \left ( \frac{\partial u_i}{\partial x_j} +\frac{\partial u_j}{\partial x_i}\right )  in \sigma_{ij}=-p\delta_{ij}+\mu\left ( \frac{\partial u_i}{\partial x_j} +\frac{\partial u_j}{\partial x_i}\right ), which is in Cartesian coordinates, to spherical coordinates. Let i =1 and j = 2,

\tau_{12}=\mu\left ( \frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1} \right )\; \; \; \; \; \; \; (53b)

In Cartesian coordinates, \nabla=\hat{x_1}\frac{\partial }{\partial x_1}+\hat{x_2}\frac{\partial }{\partial x_2}+\hat{x_3}\frac{\partial }{\partial x_3} and \textbf{\textit{u}}=u_1\hat{x_1}+u_2\hat{x_2}+u_3\hat{x_3}. So,

\left ( \frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1} \right )=\left ( \hat{x_2}\cdot \nabla \right )\left ( \textbf{\textit{u}}\cdot \hat{x_1}\right )+\left ( \hat{x_1}\cdot \nabla \right )\left ( \textbf{\textit{u}}\cdot \hat{x_2}\right )

Since \left ( \hat{x_2}\cdot \nabla \right ) and \left ( \hat{x_1}\cdot \nabla \right ) are scalars, we can rewrite the above as:

\left ( \frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1} \right )=\left [ \left ( \hat{x_2}\cdot \nabla \right )\textbf{\textit{u}}\right]\cdot \hat{x_1}+\left [ \left ( \hat{x_1}\cdot \nabla \right )\textbf{\textit{u}}\right ] \cdot \hat{x_2}\; \; \; \; \; \; \; (53c)

We can now replace the orthogonal Cartesian basis vectors with orthogonal spherical basis vectors in eq53c:

\left ( \frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1} \right )=\left [ \left ( \hat{\theta}\cdot \nabla \right )\textbf{\textit{u}}\right]\hat{r}+\left [ \left ( \hat{r}\cdot \nabla \right )\textbf{\textit{u}}\right ] \hat{\theta}

In spherical coordinates, \nabla=\frac{\partial }{\partial r}\hat{r}+\frac{1}{r}\frac{\partial }{\partial \theta}\hat{\theta}+\frac{1}{rsin\theta}\frac{\partial }{\partial \phi}\hat{\phi} and \textbf{\textit{u}}=u_r\hat{r}+u_\theta\hat{\theta}+u_\phi \hat{\phi}. So, working out the algebra and using the identities \frac{\partial \hat{r}}{\partial \theta}=\hat{\theta}, \frac{\partial \hat{\theta}}{\partial \theta}=-\hat{r} and \frac{\partial \hat{\phi}}{\partial \theta}=0, the above equation becomes:

\left ( \frac{\partial u_1}{\partial x_2}+\frac{\partial u_2}{\partial x_1} \right )=r\frac{\partial }{\partial r}\left ( \frac{u_\theta}{r} \right )+\frac{1}{r}\frac{\partial u_r}{\partial \theta}\; \; \; \; \; \; \; (53d)

Substituting eq53d in eq53b completes the exercise.

 

Substitute eq45, eq51 and r = a (on the surface of the sphere) in eq53,

\tau_{rr}=-p_\infty +\frac{3\mu\textbf{\textit{u}}}{2a}cos\theta\; \; \; \; \; \; \; (54)

Substitute eq45, eq46 and r = a (on the surface of the sphere) in eq53a,

\tau_{r\theta}=-\frac{3\mu\textbf{\textit{u}}}{2a}sin\theta\; \; \; \; \; \; \; (55)

Substitute eq54 and eq55 in eq52,

\tau=-p_\infty cos\theta+\frac{3\mu\textbf{\textit{u}}}{2a}\; \; \; \; \; \; \; (56)

The drag on the sphere FD is therefore the integral of the stress vector over the surface area of the sphere:

F_D=\int_{0}^{2\pi}\int_{0}^{\pi}\tau a^2sin\theta d\theta d\phi

Substitute eq56 in the above integral,

F_D=-2\pi a^2\int_{0}^{\pi}p_\infty sin\theta cos\theta d\theta+3\pi a\mu \textbf{\textit{u}}\int_{0}^{\pi}sin\theta d\theta

Integrate by substituting x = sinθ, dx = cosθdθ, we have

F_D=6\pi a\mu\textbf{\textit{u}}\; \; \; \; \; \; \; (57)

where a is the radius of the sphere.

 

Previous article: solution of the differential equation, E2(E2ψ) = 0
Content page of Stokes’ law
Content page of advanced chemistry
Main content page
Posted on

Continuity equation

A continuity equation in fluid mechanics is an equation that makes use of the conservation-of-mass theory to describe the transport of fluid.

With reference to the above diagram, a control volume in spherical coordinates is given by:

dV=r^2sin\theta drd\theta d\phi\;\;\;\;\;\;\;(1)

where

r is the radius of the sphere

θ is the polar angle

Φ is the azimuthal angle

Substitute the formula for density V = m/ρ in eq1 and divide throughout by ∂t,

\frac{\partial m}{\partial t}=\frac{\partial \rho}{\partial t}r^2sin\theta drd\theta d\phi\; \; \; \; \; \; \; (2)

The flow velocity of the fluid is:

\textbf{\textit{u}}=\textit u_r\hat{r}+\textit u_\theta\hat{\theta}+\textit u_\phi\hat{\phi}\; \; \; \; \; \; \; (3)

The net flow of fluid through the surfaces of the control volume in the radial direction can be expressed as the rate of change of mass through the surfaces in the radial direction:

\frac{\partial m_{out,r}}{\partial t}-\frac{\partial m_{in,r}}{\partial t}=(\rho u_r+\frac{\partial \rho u_r}{\partial r}dr)A_{out,r}-\rho u_rA_{in,r}\; \; \; \; \; \; \; (4)

Aout,r is the outflow area of the control volume in the radial direction, given by:

A_{out,r}=average\:length_{out}\:\times \:breadth_{out}

A_{out,r}=\frac{(r+dr)sin(\theta +d\theta)d\phi+(r+dr)sin\theta d\phi}{2}(r+dr)d\theta\; \; \; \; \; \; \; (5)

Substituting  sin(θ+dθ) ≈ sinθ + cosθdθ (from the Maclaurin series of cosθ and sinθ for small θ) in eq5 and eliminating higher order terms (i.e. terms with dr2, dθ2):

A_{out,r}=r^2sin\theta d\theta d\phi+2rsin\theta drd\theta d\phi\; \; \; \; \; \; \; (6)

Similarly, the inflow area of the control volume in the radial direction Ain,r is:

A_{in,r}=r^2sin\theta d\theta d\phi\; \; \; \; \; \; \; (7)

Substituting eq6 and eq7 in eq4 and eliminating higher order terms,

\frac{\partial m_{out,r}}{dt}-\frac{\partial m_{in,r}}{dt}=2\rho u_rrsin\theta drd\theta d\phi + r^2sin\theta \frac{\partial \rho u_r}{dr}drd\theta d\phi\; \; \; \; \; \; \; (8)

Repeating the above steps for the net flow of fluid through the surfaces of the control volume in the polar and azimuthal directions, we have:

\frac{\partial m_{out,\theta}}{dt}-\frac{\partial m_{in,\theta}}{dt}=\rho u_\theta rcos\theta drd\theta d\phi+ rsin\theta \frac{\partial \rho u_\theta}{\partial \theta}drd\theta d\phi \; \; \; \; \; \; \; (9)

\frac{\partial m_{out,\phi}}{dt}-\frac{\partial m_{in,\phi}}{dt}=r \frac{\partial \rho u_\phi}{\partial \phi}drd\theta d\phi \; \; \; \; \; \; \; (10)

Assuming steady fluid flow and the conservation of fluid mass, the fluxes of mass into the control volume must be equal to the fluxes of mass out of the control volume and the accumulation of mass in the control volume:

\frac{\partial m_{in,r}}{\partial t}+\frac{\partial m_{in,\theta}}{\partial t}+\frac{\partial m_{in,\phi}}{\partial t}=\frac{\partial m_{out,r}}{\partial t}+\frac{\partial m_{out,\theta}}{\partial t}+\frac{\partial m_{out,\phi}}{\partial t}+\frac{\partial m}{\partial t}\; \; \; \; \; \; \; (11)

Substituting eq8, eq9 and eq10 in eq11, we have the continuity equation in spherical coordinates.

\frac{\partial \rho}{\partial t}+\frac{1}{r^2}\frac{\partial \rho r^2u_r}{\partial r}+\frac{1}{rsin\theta} \frac{\partial \rho u_\theta sin\theta}{\partial \theta}+\frac{1}{rsin\theta} \frac{\partial \rho u_\phi}{\partial \phi}=0\; \; \; \; \; \; \; (12)

where

\frac{1}{rsin\theta} \frac{\partial \rho u_\theta sin\theta}{\partial \theta}=\frac{1}{r}\frac{\partial \rho u_\theta}{\partial \theta}+\rho u_\theta \frac{cos\theta}{rsin\theta}

 

 

\frac{1}{r^2} \frac{\partial \rho r^2u_r}{\partial r}=\frac{2}{r}\rho u_r+\frac{\partial \rho u_r}{\partial r}

 

 

NExt article: Stokes stream function
Previous article: Overview
Content page of Stokes’ law
Content page of advanced chemistry
Main content page

 

Posted on

Derivation of the differential equation: \(E^2(E^2\psi)=0\)

The derivation of the differential equation \(E^2(E^2\psi)=0\) is part of the steps in deriving Stoke’s law.

We have, in the previous section, derived the Navier-Stokes equation:

\rho\frac{D\textbf{\textit{u}}}{Dt}=-\nabla p+\rho g+\mu \nabla^2\textbf{\textit{u}}\; \; \; \; \; \; \; (32)

The terms \rho \frac{D\textbf{\textit{u}}}{Dt} and \rho g are attributed to inertia forces, while the terms \nabla p and \mu \nabla^2\textbf{\textit{u}} are due to hydrostatic forces and viscous forces respectively. For a very viscous fluid, the viscosity of the fluid is very high but fluid velocities are very low. Therefore, the non-inertia forces dominate the inertia forces and eq37 is reduced to:

\mu \nabla^2\textbf{\textit{u}}=\nabla p\; \; \; \; \; \; \; (33)

 

 

Taking the curl on both sides of eq33,

\nabla \times \mu \nabla^2\textbf{\textit{u}}=\nabla \times \nabla p

Since the curl of the gradient of a function is zero and the viscosity of an incompressible fluid is constant,

\nabla\times\mu \nabla^2\textbf{\textit{u}}=0\; \; \; \; \Rightarrow \; \; \; \; \nabla\times\nabla^2\textbf{\textit{u}}=0\; \; \; \; \; \; \; (33a)

Substituting eq14 where \nabla\cdot\textbf{\textit{u}}=0 in the vector identity \nabla^2\textbf{\textit{u}}=\nabla\left ( \nabla\cdot \textbf{\textit{u}}\right )-\nabla \times \left ( \nabla \times \textbf{\textit{u}}\right ), we get:

\nabla^2\textbf{\textit{u}}=-\nabla\times\left ( \nabla\times \textbf{\textit{u}}\right )\; \; \; \; \ (33b)

Substitute eq33b in eq33a

\nabla \times \nabla \times \left ( \nabla \times \textbf{\textit{u}}\right )=0\; \; \; \; \; \; \; (34)

Assuming that the fluid is flowing in the e_\phi direction, eq3 becomes \textbf{\textit{u}}=u_r\hat{r}+u_\theta \hat{\theta}. Noting that the curl of a function in spherical coordinates is:

\nabla\times\textbf{\textit{u}}=\frac{1}{rsin\theta}\left [ \frac{\partial }{\partial \theta}\left ( u_\phi sin\theta \right )-\frac{\partial }{\partial \phi}u_\theta \right ]\hat{r}

+\frac{1}{r}\left [\frac{1}{sin\theta} \frac{\partial }{\partial \phi} u_r-\frac{\partial }{\partial r}\left (ru_\phi \right ) \right ]\hat{\theta}+\frac{1}{r}\left [ \frac{\partial }{\partial r}\left (r u_\theta \right )-\frac{\partial }{\partial \theta}u_r \right ]\hat{\phi}

and substituting eq15 and eq16 in \textbf{\textit{u}}=u_r\hat{r}+u_\theta \hat{\theta}, we get:

\textbf{\textit{u}}=\frac{1}{r^2sin\theta}\frac{\partial \psi}{\partial \theta}\hat{r}-\frac{1}{rsin\theta}\frac{\partial \psi}{\partial r}\hat{\theta}=\nabla\times \frac{\psi}{rsin\theta}\hat{\phi}\; \; \; \; \; \; \; (35)

 

Question

Show that \nabla\times\frac{\psi}{rsin\theta}\hat{\phi}=\frac{1}{r^2sin\theta}\frac{\partial \psi}{\partial \theta}\hat{r}-\frac{1}{rsin\theta}\frac{\partial \psi}{\partial r}\hat{\theta}.

Answer

Using the curl of a function in spherical coordinates,

\nabla\times\left ( 0\hat{r}+0\hat{\theta}+\frac{\psi}{rsin\theta}\hat{\phi} \right )=\frac{1}{rsin\theta}\left [ \frac{\partial }{\partial \theta}\left ( \frac{\psi}{rsin\theta} sin\theta\right )-\frac{\partial }{\partial \phi}0 \right ]\hat{r}

+\frac{1}{r}\left [ \frac{1}{sin\theta}\frac{\partial }{\partial \phi}0-\frac{\partial }{\partial r}\left ( \frac{r\psi}{rsin\theta} \right ) \right ]\hat{\theta}\: +\frac{1}{r}\left [ \frac{\partial }{\partial r}r0-\frac{\partial }{\partial \theta}0 \right ]\hat{\phi}\:

\nabla\times\frac{\psi}{rsin\theta}\hat{\phi}=\frac{1}{r^2sin\theta}\frac{\partial \psi}{\partial \theta}\hat{r}-\frac{1}{rsin\theta}\frac{\partial \psi}{\partial r}\hat{\theta}

 

Substituting eq35 in eq34 and working out the algebra, we have:

\nabla\times\nabla\times\left ( -\frac{1}{rsin\theta}E^2\psi\hat{\phi} \right )=0\; \; \; \; \; \; \; (35a)

where E^2=\frac{\partial ^2}{\partial r^2}+\frac{sin\theta}{r^2}\frac{\partial }{\partial \theta}\left ( \frac{1}{sin\theta}\frac{\partial }{\partial \theta} \right ).

Comparing eq35a with eq34,

\nabla\times\textbf{\textit{u}}=-\frac{1}{rsin\theta}E^2\psi\hat{\phi}\; \; \; \; \; \; \; (35b)

Continuing with the algebra for eq35a, we get:

E^2\left ( E^2\psi \right )=0\; \; \; \; \; \; \; (36)

 

NExt article: solution of differential equation
Previous article: navier-stokes equation
Content page of Stokes’ law
Content page of advanced chemistry
Main content page
Posted on

Navier-Stokes equation (derivation)

The Navier-Stokes equation is used to describe the motion of a viscous fluid and was named after Claude-Louis Navier, a French scientist, and George Stokes, an Irish scientist. To derive the equation, let’s consider the flow of an incompressible fluid through an infinitesimal cubic control volume* depicted in the diagram below.

* the control volume here refers to the volume of an object in the fluid.

The x-component of the velocity of the fluid is a function of x, y, z and t, i.e. ux (x, y, z, t) where t is time. The total differential of ux is:

du_x=\frac{\partial u_x}{\partial x}dx+\frac{\partial u_x}{\partial y}dy+\frac{\partial u_x}{\partial z}dz+\frac{\partial u_x}{\partial t}dt

Dividing the equation by dt and letting dx/dt = ux, dy/dt = uy, dz/dt = uz, we have:

\frac{du_x}{dt}=u_x\frac{\partial u_x}{\partial x}+u_y\frac{\partial u_x}{\partial y}+u_z\frac{\partial u_x}{\partial z}+\frac{\partial u_x}{\partial t}\; \; \; \; \; \; \; (24)

Let’s analyse the hydrostatic and hydrodynamic forces acting on each face of the infinitesimal control volume in the x-direction (see diagram above where point C is the origin). They are summarised as follows:

Face

 Force per unit area

Force

BDFH

 σxx(x + dx) σxx(x + dx)dydz

ACEG

 σxx(x) σxx(x)dydz

ABHG

 τyx(y + dy) τyx(y + dy)dxdz

CDFE

 τyx(y) τyx(y)dxdz

ABDC

 τzx(z) τzx(z)dxdy

GHFE

 τzx(z + dz) τzx(z + dz)dxdy

σxx = F/A , i.e. the force F acting on a unit area of the x-plane in the x-direction, thereby the double-x subscript. The notation σxx(x) refers to σxx evaluated at a distance of (x) while σxx(x + dx) is evaluated at a distance of (x + dx). A different notation τ is given to the forces that act parallel to the faces of the control volume, e.g. τyx(y + dy) is the force acting on a unit area of the y-plane in the x-direction that is evaluated at a distance of (y + dy).

With reference to the defined directions of the forces, FH the sum of hydrostatic and hydrodynamic forces acting on the infinitesimal control volume in the positive x-direction is:

In addition to hydrostatic and hydrodynamic forces, an inertia force Fg (gravitational force) also acts on the infinitesimal control volume in the x-direction:

F_g=\rho g_xdxdydz

where ρ is the density of the fluid and gx is the x-component of the acceleration due to gravity.

The total force FT acting on the infinitesimal control volume in the x-direction is therefore:

F_T=F_H+F_g\; \; \; \; \; \; \; (25)

Next, by i) substituting FH and Fg and eq24 in eq25; ii) dividing the substituted equation by dxdydz and iii) letting dx, dy, dz → 0 for an infinitesimal control volume, we have:

\rho\left ( u_x\frac{\partial u_x}{\partial x} +u_y\frac{\partial u_x}{\partial y}+u_z\frac{\partial u_x}{\partial z} +\frac{\partial u_x}{\partial t} \right )=\frac{\partial \sigma_{xx}}{\partial x} +\frac{\partial \tau_{yx}}{\partial y} +\frac{\partial \tau_{zx}}{\partial z} +\rho g_x\; \; \; \; \; \; \; (26)

Using the same logic and repeating the above steps, we have for the y-direction and the z-direction:

\rho\left ( u_x\frac{\partial u_x}{\partial x}+ u_y\frac{\partial u_x}{\partial y}+ u_z\frac{\partial u_x}{\partial z}+\frac{\partial u_y}{\partial t}\right )= \frac{\partial \sigma_{yy}}{\partial y}+\frac{\partial \tau_{zy}}{\partial z}+\frac{\partial \tau_{xy}}{\partial x}+\rho g_y\; \; \; \; \; \; \; (27)

\rho\left ( u_x\frac{\partial u_x}{\partial x}+ u_y\frac{\partial u_x}{\partial y}+ u_z\frac{\partial u_x}{\partial z}+\frac{\partial u_z}{\partial t}\right )= \frac{\partial \sigma_{zz}}{\partial z}+\frac{\partial \tau_{xz}}{\partial x}+\frac{\partial \tau_{yz}}{\partial y}+\rho g_z\; \; \; \; \; \; \; (28)

From the articles on constitutive relation, \sigma_{ij}=-p\delta_{ij}+\mu\left ( \frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i}\right ). We therefore substitute \sigma_{xx}=-p+2\mu \frac{\partial u_x}{\partial x} , \tau_{yx}=\mu\left ( \frac{\partial u_x}{\partial y}+\frac{\partial u_y}{\partial x}\right ) and \tau_{zx}=\mu\left ( \frac{\partial u_z}{\partial x}+\frac{\partial u_x}{\partial z}\right ) in eq26, noting that \frac{\partial u_x}{\partial x}+\frac{\partial u_y}{\partial y}+\frac{\partial u_z}{\partial z}=\nabla\cdot\textbf{\textit{u}}=0 (see eq14) to give:

\rho\left ( \frac{\partial u_x}{\partial t}+u_x\frac{\partial u_x}{\partial x} +u_y\frac{\partial u_x}{\partial y}+u_z\frac{\partial u_x}{\partial z} \right )=\rho g_x-\frac{\partial p}{\partial x}+\mu\left ( \frac{\partial^2 u_x}{\partial x^2} +\frac{\partial^2 u_x}{\partial y^2}+\frac{\partial^2 u_x}{\partial z^2}\right )\; \; \; \; \; \; \; (29)

Using the same logic and repeating the above step for eq27 and eq28, we have:

\rho\left ( \frac{\partial u_y}{\partial t} +u_x\frac{\partial u_y}{\partial x} +u_y\frac{\partial u_y}{\partial y}+u_z\frac{\partial u_y}{\partial z} \right )=\rho g_y-\frac{\partial p}{\partial y}+\mu\left ( \frac{\partial ^2u_y}{\partial x^2}+\frac{\partial ^2u_y}{\partial y^2} +\frac{\partial ^2u_y}{\partial z^2}\right )\; \; \; \; \; \; \; (30)

\rho\left ( \frac{\partial u_z}{\partial t} +u_x\frac{\partial u_z}{\partial x} +u_y\frac{\partial u_z}{\partial y}+u_z\frac{\partial u_z}{\partial z} \right )=\rho g_z-\frac{\partial p}{\partial z}+\mu\left ( \frac{\partial ^2u_z}{\partial x^2}+\frac{\partial ^2u_z}{\partial y^2} +\frac{\partial ^2u_z}{\partial z^2}\right )\; \; \; \; \; \; \; (31)

Eq29, eq30 and eq31 are the Cartesian form of the Navier-Stokes equations, which when multiplied throughout by the unit vectors \hat{x}, \hat{y} and \hat{z} respectively and summed, becomes the vector form:

\rho\frac{D\textbf{\textit{u}}}{Dt}=-\nabla p+\rho g+\mu \nabla^2\textbf{\textit{u}}\; \; \; \; \; \; \; (32)

where

\nabla^2\textbf{\textit{u}}=\nabla^2 u_x\hat{x}+\nabla^2 u_y\hat{y}+\nabla^2 u_z\hat{z}

\frac{D}{Dt} is the convective derivative defined as \frac{D}{Dt}=\frac{\partial }{\partial t}+\textbf{\textit{u}}\cdot \nabla

\nabla is the gradient defined as \nabla =\hat{x}\frac{\partial }{\partial x}+\hat{y}\frac{\partial }{\partial y}+\hat{z}\frac{\partial }{\partial z}

\nabla^2 is the vector Laplacian defined as \nabla^2\textbf{\textit{u}}=\nabla^2 u_x\hat{x}+\nabla^2 u_y\hat{y}+\nabla^2 u_z\hat{z}

\mu is the viscosity of the fluid

 

NExt article: Derivation of differential equation, E2(E2ψ) = 0
Previous article: stokes stream function
Content page of Stokes’ law
Content page of advanced chemistry
Main content page
Posted on

Solution of the differential equation: \(E^2(E^2\psi)=0\)

How do we solve the differential equation \(E^2(E^2\psi)=0\)?

We have, in an earlier article, defined the Stokes stream function for an incompressible fluid at the boundaries of r = a and r = ∞ where:

\psi=0\; \; \; \; \; if\; r=a

\psi=\frac{\textbf{\textit{u}}}{2}r^2sin^2\theta\; \; \; \; \; if\; r=\infty

To develop a Stokes stream function that satisfies eq36 for a ≤ r ≤ ∞ , let’s consider a solution of the form:

\psi=fsin^2\theta\; \; \; \; \; \; \; (37)

where f is a function of r.

Substitute eq37 in eq36 and noting that \frac{\partial }{\partial \theta}\left ( \frac{1}{sin\theta}\frac{\partial }{\partial \theta} \right )sin^2\theta=-2sin\theta, we have:

\left ( \frac{\partial^2 }{\partial r^2} -\frac{2}{r^2}\right )^2f=0\; \; \; \; \; \; \; (38)

Letting f=r^n and solving eq38, we get n = ±1, 2, 4. Therefore, the functionftakes the form:

f=\frac{A}{r}+Br+Cr^2+Dr^4\; \; \; \; \; \; \; (39)

Substitute eq39 in eq37,

\psi=\left ( \frac{A}{r}+Br+Cr^2+Dr^4 \right )sin^2\theta\; \; \; \; \; \; \; (40)

To satisfy the function at the top of the page at r → ∞, we let D = 0 and C = u/2. Eq40 becomes:

\psi=\left ( \frac{A}{r}+Br+\frac{\textbf{\textit{u}}}{2}r^2 \right )sin^2\theta\; \; \; \; \; \; \; (41)

Substitute eq41 in eq15 and eq16, we have:

u_r=2cos\theta\left ( \frac{A}{r^3}+\frac{B}{r}+\frac{\textbf{\textit{u}}}{2} \right )\; \; \; \; \; \; \; (42)

u_\theta=sin\theta\left (- \frac{A}{r^3}+\frac{B}{r}+\textbf{\textit{u}} \right )\; \; \; \; \; \; \; (43)

To satisfy the function at the top of the page at r = a, ur = 0 and uθ = 0, i.e. velocity of the fluid on the surface of the object is zero relative to the velocity of the object. This is because the force of attraction between the fluid particles and solid particles of the object is greater than that between the fluid particles. Eq42 and eq43 become:

\frac{A}{a^3}+\frac{B}{a}+\frac{\textbf{\textit{u}}}{2}=0\; \; \; \; and\; \; \; \;-\frac{A}{a^3}+\frac{B}{a}+\textbf{\textit{u}}=0

Solving the simultaneous equations gives:

A=\frac{\textbf{\textit{u}}}{4}a^3\; \; \; \; and\; \; \; \; B=-\frac{3\textbf{\textit{u}}}{4}a

Substitute A and B in eq41,

\psi=sin^2\theta\frac{\textbf{\textit{u}}}{4}\left ( \frac{a^3}{r}-3ar+2r^2 \right )\; \; \; \; \; \; \; (44)

Therefore, we have a Stokes stream function that satisfies the differential equation eq36 for a ≤ r ≤ ∞ .

Substitute eq44 in eq15 and eq16, we have:

u_r=\textbf{\textit{u}}cos\theta\left ( \frac{a^3}{2r^3}-\frac{3a}{2r}+1 \right )\; \; \; \; \; \; \; (45)

u_\theta=-\textbf{\textit{u}}sin\theta\left ( -\frac{a^3}{4r^3}-\frac{3a}{4r}+1 \right )\; \; \; \; \; \; \; (46)

Eq45 and eq46 are the component flow velocities of the viscous incompressible fluid.

 

NExt article: derivation of stokes’ law
Previous article: derivation of the differential equation, E2(E2ψ) = 0
Content page of Stokes’ law
Content page of advanced chemistry
Main content page
Posted on

Stokes stream function

The Stokes stream function is a mathematical representation of the trajectories of particles in a steady flow of fluid over an object. In other words, a plot of the Stokes stream function results in streamlines seen in the diagram below.

To derive the Stokes stream function, we make the following assumptions:

  • The flow of the fluid is axisymmetric, i.e. defined along the polar axis (z-axis), and therefore the velocity components are independent of Φ. 
  • The fluid is incompressible and therefore has a constant density.

Eq12 of the previous article becomes:

\frac{1}{r^2}\frac{\partial r^2u_r}{\partial r}+\frac{1}{rsin\theta}\frac{\partial u_\theta sin\theta}{\partial \theta}=0\; \; \; \; \; \; \; (13)

 

Question

Show that \nabla\cdot u=0.

Answer

The divergence of a function in spherical coordinate is:

\nabla\cdot \boldsymbol{\mathit{f}}=\frac{1}{r^{2}}\frac{\partial r^{2}f_r}{\partial r}+\frac{1}{rsin\theta}\frac{\partial f_{\theta}sin\theta}{\partial\theta}+\frac{1}{rsin\theta}\frac{\partial f_{\phi}}{\partial \phi}

If the flow is axisymmetric, the divergence of an axisymmetric flow function in spherical coordinates is exactly the LHS of eq13 and we can write eq13 as:

\nabla \cdot\textbf{\textit{u}}=0\; \; \; \; \; \; \; (14)

Eq14 is needed later for the derivation of the differential equation E2(E2ψ=0).

 

George Stokes developed the solution to eq13 by defining the Stokes stream function ψ where:

u_r=\frac{1}{r^2sin\theta}\frac{\partial \psi}{\partial \theta}\; \; \; \; \; \; \; (15)

and

u_\theta=-\frac{1}{rsin\theta}\frac{\partial \psi}{\partial r}\; \; \; \; \; \; \; (16)

 

Question

Show that eq15 and eq16 satisfy eq13.

Answer

Substitute eq15 and eq16 in eq13:

\frac{1}{r^2sin\theta}\frac{\partial }{\partial r}\left ( \frac{\partial \psi}{\partial \theta} \right )-\frac{1}{r^2sin\theta}\frac{\partial }{\partial \theta}\left ( \frac{\partial \psi}{\partial r} \right )=0

Since \frac{\partial }{\partial r}\left ( \frac{\partial \psi}{\partial \theta} \right )= \frac{\partial }{\partial \theta}\left ( \frac{\partial \psi}{\partial r} \right ) , eq15 and eq16 satisfy eq13.

 

With reference to the diagram at the top of the page, the flow velocity u, which is defined in the z-direction, varies at different distances from the surface of the sphere. On the surface of the sphere,

r=a\; \; \; \; and\; \; \; \; \psi=0\; \; \; \; \; \; \; (17)

At r = ∞, we assume that the flow velocity of the fluid is uniform. The flow velocity at any point in fluid at can be deconstructed into its radial and polar components (see above diagram) where:

u_r=\textbf{\textit{u}}cos\theta\; \; \; \; \; \; \; (18)

and

u_\theta=-\textbf{\textit{u}}sin\theta\; \; \; \; \; \; \; (19)

Substitute eq15 in eq18 and eq16 in eq19, we have,

\frac{1}{r^2sin\theta}\frac{\partial \psi}{\partial \theta}= \textbf{\textit{u}}cos\theta\; \; \; \; \; \; \; (20)

\frac{1}{rsin\theta}\frac{\partial \psi}{\partial r}= \textbf{\textit{u}}sin\theta\; \; \; \; \; \; \; (21)

Integrating eq21, we get:

\textbf{\textit{u}}sin^2\theta\int_{a}^{r}r\; dr=\int_{0}^{\psi}d\psi

\psi=\frac{\textbf{\textit{u}}}{2}\left ( r^2-a^2 \right )sin^2\theta

At r = ∞, r2 » a2, so

\psi\approx\frac{\textbf{\textit{u}}}{2}r^2sin^2\theta\; \; \; \; \; \; \; (22)

Eq17 and eq22 express the Stokes stream function for an incompressible fluid at the boundaries of r = a and r = ∞ respectively. In the next few articles, we shall define the Stokes stream function at a ≤ r ≤ ∞.

 

NExt article: Navier-Stokes equation
Previous article: continuity equation
Content page of Stokes’ law
Content page of advanced chemistry
Main content page
Mono Quiz