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

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

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 articleE^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 Fis 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

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(θ+) ≈ 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, 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

 

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

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 uis:

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/dtuxdy/dtuydz/dtuz, 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, Fthe 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 Facting 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

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 ar ≤ ∞ , 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 ar ≤ ∞ .

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

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

Stokes’ law: overview

Stokes’ law, which was derived in 1851 by an Irish scientist named George Stokes, describes the drag force exerted on a spherical object in a viscous fluid through the formula:

F_d=6\pi\mu ru

where

Fd is the frictional force (drag) acting on the surface of the object

μ is the viscosity of the fluid

r is the radius of the spherical object

u is the flow velocity of the fluid

 

 

The derivation of the equation for Stokes’ law involves the following steps:

Step 1: Derive the continuity equation in spherical coordinates for an incompressible fluid.

Step 2:  Define the Stokes stream function to satisfy the continuity equation for an incompressible fluid.

Step 3:  Derive the Navier-Stokes equations.

Step 4:  Derive the differential equation, E2(E2ψ) = 0.

Step 5:  Solve the differential equation.

Step 6:  Derive the Stokes’ law equation.

 

next article: Continuity equation
Content page of Stokes’ law
Content page of advanced chemistry
Main content page

Reciprocal lattice

A reciprocal lattice is the diffraction image of a real crystal lattice. It is formed by linear combinations of basis vectors of the reciprocal space.

In an earlier article, we mentioned that the three Laue equations must be simultaneously satisfied for constructive interference to occur in three dimensions. In other words, the solution to the Laue equations is a mathematical description of the diffraction pattern of an X-ray diffraction experiment.

Let h = (s – s0) and we can rewrite eq20, eq21 and eq22 as:

\textbf{\textit{a}}\cdot\textbf{\textit{h}}=h\; \; \; \; \; \; \; (24a)

\textbf{\textit{b}}\cdot\textbf{\textit{h}}=k\; \; \; \; \; \; \; (24b)

\textbf{\textit{c}}\cdot\textbf{\textit{h}}=l\; \; \; \; \; \; \; (24c)

A general solution to the simultaneous equations of eq24a, eq24b and eq24c is:

\textbf{\textit{h}}=h\textbf{\textit{a}}^*+k\textbf{\textit{b}}^*+l\textbf{\textit{c}}^*\; \; \; \; \; \; \; (24d)

where \textbf{\textit{a}}^*=\frac{\textbf{\textit{b}}\times\textbf{\textit{c}}}{V}, \textbf{\textit{b}}^*=\frac{\textbf{\textit{c}}\times\textbf{\textit{a}}}{V} and \textbf{\textit{c}}^*=\frac{\textbf{\textit{a}}\times\textbf{\textit{b}}}{V}. V is the volume of the unit cell formed by the basis vectors a, b and c.

 

Question

Show that eq24d is a general solution to the Laue equations.

Answer

Substitute eq24d in eq24a

\textbf{\textit{a}} \cdot (h\textbf{\textit{a}}^*+k\textbf{\textit{b}}^*+l\textbf{\textit{c}}^*)=h

h\frac{\textbf{\textit{a}}\cdot (\textbf{\textit{b}}\times \textbf{\textit{c}})}{V} +k\frac{\textbf{\textit{a}}\cdot (\textbf{\textit{c}}\times \textbf{\textit{a}})}{V}+l\frac{\textbf{\textit{a}}\cdot (\textbf{\textit{a}}\times \textbf{\textit{b}})}{V}=h

Since (c × a) is a vector that is perpendicular to the plane formed by c and a, a · (c × a) = IaI Ic × aI cos90º = 0. By the same logic, a · (a × b) = 0. Therefore,

h\frac{\textbf{\textit{a}}\cdot (\textbf{\textit{b}}\times \textbf{\textit{c}})}{V} =h

Since a · (b × c) = V (see below for proof), LHS = RHS. We can substitute eq24d in eq24a and eq24b with similar results.

 

The vector h in eq24d is also a lattice vector but with basis vectors of a*, b* and c* instead of a, b and c. The dimension of a* is length-1 because \textbf{\textit{a}}^*=\frac{\textbf{\textit{b}}\times\textbf{\textit{c}}}{V}=\frac{ \left |\textbf{\textit{b}} \right | \left | \textbf{\textit{c}}\right |sin\theta\textbf{\textit{n}}}{V}, where n is a unit vector in the direction of a(the unit vector n is necessary because a* is a vector but \frac{\left | \textbf{\textit{b}}\right |\left | \textbf{\textit{c}}\right |sin\theta }{V} is a scalar). Similarly, b* and c* have dimensions of the reciprocal of length. Hence, h is known as a reciprocal lattice vector.

Just as a lattice vector describes the positions of lattice points in a space lattice, the reciprocal lattice vector mathematically defines the positions of lattice points in a reciprocal space lattice, which means that the diffraction pattern is constructed in a reciprocal space (see diagram below).

Even though the reciprocal lattice and the Ewald sphere allow crystallographers to visualise diffraction patterns in relation to various X-ray diffraction techniques, Bragg’s law is most often used to explain X-ray diffraction concepts. This is because Bragg’s law is easier to apply than the Laue equations, and in some cases involves the knowledge of just an angle θ to deduce the unit cell types and unit cell dimensions of a sample.

 

Question

Show that a · (b × c) = V.

Answer

With reference to the above diagram, a, b, c are the basis vectors of the unit cell. The area of the base of the unit cell, which is a parallelogram, is:

base\: area=\vert\boldsymbol{\mathit{b}}\vert\vert\boldsymbol{\mathit{c}}\vert sin\theta=\vert\boldsymbol{\mathit{b}}\times\boldsymbol{\mathit{c}}\vert

The volume of the unit cell, V, is:

V=(base\: area)(height)=\left | \textbf{\textit{b}}\times\textbf{\textit{c}}\right | \left |\textbf{\textit{a}} \right |cos\phi=\textbf{\textit{a}}\cdot (\textbf{\textit{b}}\times\textbf{\textit{c}})

a · (b × c) is known as the scalar triple product.

 

next article: ewald sphere
previous article: equivalence of the laue equations and the bragg equation
Content page of X-ray crystallography
Content page of advanced chemistry
Main content page
Mono Quiz