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

Electron density as a Fourier transform of the structure factor

The electron densities of atoms in a crystal are Fourier transforms of the structure factor.

 

In the article on scattering factor, we have restricted the electron density ρ to a lattice point, which is too simplistic. We have also described the scattering factor (eq27) as

df=\rho e^{i\phi}dV

or in its integrated form

f=\int \rho e^{i\phi}dV\; \; \; \; \; \; \; (40)

If we now extend the distribution of ρ through the unit cell, eq40 becomes the structure factor:

F_{hkl}=\int \rho 'e^{i\phi}dV

where ρ’ρ(xyz) is the electron density at coordinates xyz in the unit cell and \phi=2\pi(\frac{hx}{a}+\frac{ky}{b}+\frac{lz}{c} ), i.e.

F_{hkl}=\int \rho(xyz)e^{i2\pi(\frac{hx}{a}+\frac{ky}{b}+\frac{lz}{c} )}dV\; \; \; \; \; \; \; (41)

Let dV be an infinitesimal volume of the unit cell with edges dx, dy, dz that are parallel to the unit cell axes of a, b, c (volume of unit cell is V). The ratio of dV/V must be equal to (dxdydz)/abc and we can rewrite eq41 as

F_{hkl}=\frac{V}{abc}\int_{0}^{a}\int_{0}^{b}\int_{0}^{c}\rho (xyz)e^{i2\pi(\frac{hx}{a}+\frac{ky}{b}+\frac{lz}{c})}dxdydz\; \; \; \; \; \; \; (42)

A Fourier series is an expansion series used to represent a periodic function and is given by:

f(x)=\sum_{n=0}^{\infty }[A_ncos(nx)+B_nsin(nx)]

or by its complex form

f(x)=\sum_{n=-\infty }^{\infty }c_ne^{i\frac{2\pi nx}{a}}

Due to the repetitive arrangement of atoms in a crystal, the electron density in a crystal is also periodic and can be expressed as a Fourier series:

\rho (x)=\sum_{n=-\infty }^{\infty }c_ne^{i\frac{2\pi nx}{a}}\; \; \; \; \; \; \; (43)

In three dimensions, eq43 becomes

\rho (xyz)=\sum_{n=-\infty }^{\infty }\sum_{m=-\infty }^{\infty }\sum_{o=-\infty }^{\infty }c_{nmo}e^{i2\pi(\frac{nx}{a}+\frac{ my}{b}+\frac{ oz}{c})}\; \; \; \; \; \; \; (44)

Substitute eq44 in eq42

F_{hkl}=\frac{V}{abc}\int_{0}^{a}\int_{0}^{b}\int_{0}^{c}\sum_{n=-\infty }^{\infty }\sum_{m=-\infty }^{\infty }\sum_{o=-\infty }^{\infty }c_{nmo}e^{i2\pi(\frac{nx}{a}+\frac{my}{b}+\frac{oz}{c})} e^{i2\pi(\frac{hx}{a}+\frac{ky}{b}+\frac{lz}{c})}dxdydz

F_{hkl}=\frac{V}{abc}\sum_{n=-\infty }^{\infty }\sum_{m=-\infty }^{\infty }\sum_{o=-\infty }^{\infty }c_{nmo}\int_{0}^{a}e^{i2\pi\frac{x}{a}(n+h)}dx\int_{0}^{b}e^{i2\pi\frac{y}{b}(m+k)}dy\int_{0}^{c}e^{i2\pi\frac{z}{c}(o+l)}dz\; \; \; \; (45)

If n ≠ –h,

\int_{0}^{a}e^{i2\pi\frac{x}{a}(n+h)}dx=\int_{0}^{a}[cos2\pi\frac{x}{a}(n+h)+isin2\pi\frac{x}{a}(n+h)]dx

=\int_{0}^{a}cos2\pi\frac{x}{a}(n+h)dx=0

which makes Fhkl = 0.

Similarly, Fhkl = 0 if m ≠ -k or o ≠ -l. Therefore, the surviving term in the triple summation in eq45 corresponds to the case of n = –hm = –k, o = –l, giving

F_{hkl}=\frac{V}{abc}c_{-h,-k,-l}\int_{0}^{a}dx\int_{0}^{b}dy\int_{0}^{c}dz=Vc_{-h,-k,-l}

c_{-h,-k,-l}=\frac{F_{hkl}}{V}\; \; \; \; \; \; \; (46)

When n = –hm = –k, o = –l, eq44 becomes

\rho (xyz)=\sum_{-h=-\infty }^{\infty }\sum_{-k=-\infty }^{\infty }\sum_{-l=-\infty }^{\infty } c_{-h,-k,-l}e^{-i2\pi(\frac{hx}{a}+\frac{ky}{b}+\frac{lz}{c})}

Since \sum_{-h=-\infty }^{\infty }\sum_{-k=-\infty }^{\infty }\sum_{-l=-\infty }^{\infty }=\sum_{h=\infty }^{-\infty }\sum_{k=\infty }^{-\infty }\sum_{l=\infty }^{-\infty }=\sum_{h=-\infty }^{\infty }\sum_{k=-\infty }^{\infty }\sum_{l=-\infty }^{\infty }

\rho (xyz)=\sum_{h=-\infty }^{\infty }\sum_{k=-\infty }^{\infty }\sum_{l=-\infty }^{\infty }c_{-h,-k,-l}e^{-i2\pi(\frac{hx}{a}+\frac{ky}{b}+\frac{lz}{c})}\; \; \; \; \; \; \; (47)

Substitute eq46 in eq47

\rho (xyz)=\frac{1}{V}\sum_{h=-\infty }^{\infty }\sum_{k=-\infty }^{\infty }\sum_{l=-\infty }^{\infty }F_{hkl}e^{-i2\pi(\frac{hx}{a}+\frac{ky}{b}+\frac{lz}{c})}\; \; \; \; \; \; \; (48)

Since the limits of the integrals in eq42 can be changed to run from to without affecting the values of the integrals because is zero outside the crystal, we have

Eq48 and eq48a are a Fourier transform pair. As mentioned in the previous article, the intensity of a diffraction signal is proportional to the square of the magnitude of the three-dimensional structure factor, i.e. I \propto \left | F_{hkl} \right |^2. If we know the value of Fhkl (which in principle is the square root of the intensity of a peak from an X-ray diffraction experiment \sqrt{\left | F_{hkl} \right |^2}=\left | F_{hkl} \right | ) and having indexed the plane contributing to this intensity peak (i.e. knowing the h, k, l values), we can determine ρ(xyz) using a mathematical software. The solution to ρ(xyz) is an electron density map that elucidates bond lengths and bond angles of the compound. However, a problem called the phase problem arises.

 

next article: phase problem
previous article: structure factor
Content page of X-ray crystallography
Content page of advanced chemistry
Main content page

Scattering factor (crystallography)

The scattering factor describes the amplitude of scattered rays from a single atom.

The presence or absence of an intensity peak depends on the interference of scattered X-rays at the detector of the powder X-ray diffractometer. Since the interference of scattered X-rays is the sum of the amplitudes of X-ray waves from many atoms in a sample, we must begin with the analysis of the resultant amplitude of scattered rays from a single atom (scattering factor), which is, in turn, the sum of amplitudes of waves scattered by the electrons in the atom.

The amplitude of a travelling wave is expressed by the equation:

y=Acos(kx-\omega t+\phi )

We can also expressed the general wave equation in its complex form:

y=Ae^{i(kx-\omega t+\phi )}\; \; \; \; \; \; \; (25)

The complex form of the wave equation is equivalent to the travelling wave equation if we consider only the real part of it: Re(y), which is the value that is physically observed in an X-ray diffraction experiment. The reason eq25 is preferred over the travelling wave equation is that it is easier to manipulate mathematically to derive the scattering factor. Recall that ω = 2π/T (angular frequency), k = 2π/λ (wave number) and Φ is the phase difference. If we assume that the scattering of X-rays by an atom is elastic, i.e. there is no change in frequency of the X-ray between the incident and scattered rays, the term kxωt becomes a constant and we can simplify eq25 to:

y=A'e^{i\phi }\; \; \; \; \; \; \; (26)

where A’ = Aei(kx-ωt).

The amplitudes of the scattered X-rays not only differ in terms of their phase as suggested by eq26 but are also proportional to the number of electrons in the atom, ρdV, where ρ and V are electron density and volume of the atom respectively. If we define A’ as the number of electrons in the atom, eq26 becomes

df=\rho e^{i\phi }dV\; \; \; \; \; \; \; (27)

where the amplitude of the scattered X-rays from an atom is denoted by f, the scattering factor of an atom.

To determine an expression for Φ, we note that the ratio of path difference δ and wavelength is equal to the ratio of phase difference and 2π, i.e.

\frac{\delta }{\lambda }=\frac{\phi }{2\pi }\; \; \; \; \; \; \; (28)

Substitute eq19 in eq28,

\phi =2\pi\textbf{\textit{a}}\cdot (\textbf{\textit{s}}-\textbf{\textit{s}}_0)\; \; \; \; \; \; \; (29)

Since \textbf{\textit{a}}\cdot (\textbf{\textit{s}}-\textbf{\textit{s}}_0)=\left |\textbf{\textit{a}} \right |\left | \textbf{\textit{s}}-\textbf{\textit{s}}_0\right |cos\mu and from eq23\left | \textbf{\textit{s}}-\textbf{\textit{s}}_0\right |=\frac{2sin\theta }{\lambda }, eq29 becomes:

\phi =\frac{4\pi }{\lambda }\left | \textbf{\textit{a}}\right |sin\theta cos\mu

In three dimensions, we let a = r and IaI = IrI = r,

\phi =\frac{4\pi }{\lambda }rsin\theta cos\mu\; \; \; \; \; \; \; (30)

Substitute eq30 in eq27

df=\rho e^{ikrcos\mu }dV\; \; \; \; \; \; \; (31)

where k=\frac{4\pi }{\lambda }sin\theta

Since ρ is spherically symmetrical, we can represent the volume element dV in eq31 in spherical coordinates where dV = r2sinμdrdμdΦ (see above diagram). Eq31 becomes:

df=\rho e^{ikrcos\mu }r^2sin\mu drd\mu d\phi

Integrating both sides,

f=\int_{0}^{\infty }\int_{0}^{\pi }\int_{0}^{2\pi }\rho e^{ikrcos\mu }r^2sin\mu drd\mu d\phi

f=2\pi \int_{0}^{\infty }\rho r^2dr\int_{0}^{\pi }e^{ikrcos\mu }sin\mu d\mu \; \; \; \; \; \; \; (32)

Note that d(cosμ) = -sinμdμ; when μ = π, cosμ = -1; when μ = 0, cosμ = 1. So the second integral in eq32 becomes

\int_{0}^{\pi }e^{ikrcos\mu }sin\mu d\mu =-\int_{1}^{-1 }e^{ikrcos\mu } d(cos\mu)=\int_{-1}^{1 }e^{ikrcos\mu }d(cos\mu)

=\left [ \frac{1}{ikr}e^{ikrcos\mu } \right ]_{-1}^{1}=\frac{e^{ikr}-e^{-ikr}}{ikr}

Eq32 now becomes

f=2\pi \int_{0}^{\infty }\rho r^2 \frac{e^{ikr}-e^{-ikr}}{ikr}dr

Since e^{ikr}-e^{-ikr}=coskr+isinkr-(coskr-isinkr)=2isinkr

f=4\pi \int_{0}^{\infty }\rho \frac{sinkr}{kr}r^2dr\; \; \; \; \; \; \; (33)

where k=\frac{4\pi }{\lambda }sin\thetaρ is a function of r and remains within the integral.

So far, we have developed an expression (eq33) that describes the resultant amplitude of scattered rays from a single atom. To further analyse the interference of scattered X-rays from multiple atoms in a sample, we have to derive another expression called the structure factor.

 

next article: structure factor
previous article: single crystal X-ray diffraction
Content page of X-ray crystallography
Content page of advanced chemistry
Main content page

Ewald sphere (crystallography)

The Ewald sphere is a mathematical construct that relates the geometry of the incident and scattered wave vectors to the reciprocal lattice.

Consider a crystal at A being irradiated by an incident wave vector s0, scattering a wave vector s that makes an angle 2θ with s0, thereby satisfying Bragg’s law (see diagram below).

Since X-ray scattering is elastic, IsI = Is0I and the two wave vectors become radii of a sphere called the Ewald sphere. The vector OP is the reciprocal lattice vector (s – s0) that is denoted by h with the origin at O.

\frac{OP}{IO}=\frac{\left | \textbf{\textit{h}}\right |}{IO}=sin\theta\; \; \; \; \; \; \; (24e)

Since h = (s – s0), eq24 becomes

d_{nh,nk,nl}=\frac{1}{\left | \textbf{\textit{h}}\right |}\; \; \; \; \; \; \; (24f)

From eq11, sin\theta =\frac{\lambda }{2d_{nh,nk,nl}}. Substitute sin\theta =\frac{\lambda }{2d_{nh,nk,nl}} and eq24f in eq24e, we have:

IO=\frac{2}{\lambda }

This means that IA = AO = AP = 1/λ. Therefore, to satisfy Bragg’s law that results in constructive interference of scattered X-rays, the head of the reciprocal lattice vector must lie on any point on the surface of the Ewald sphere. From eq24f and eq11,

\left | \textbf{\textit{h}}\right |=\frac{2sin\theta }{\lambda }

Since -1 ≤ sinθ ≤ 1,

\left | \textbf{\textit{h}}\right |_{max}=\frac{2 }{\lambda }\; \; \; \; \; \; \; (24g)

As the maximum magnitude of the reciprocal lattice vector is 2/λ, all reciprocal lattice vectors that potentially satisfy the Laue equations or Bragg’s law are enclosed in a sphere of radius known as the limiting sphere. The Ewald sphere is used to visualise different X-ray diffraction techniques including single crystal X-ray diffraction.

next article: single crystal X-ray diffraction
previous article: reciprocal lattice
Content page of X-ray crystallography
Content page of advanced chemistry
Main content page

Equivalence of the Laue equations and the Bragg equation

The Laue equations are equal to the Bragg equation if we denote the angle between the wave vectors s and s0  by 2θ (see diagram below), which makes the vector s – s0  normal to the lattice plane.

As the scattering of X-rays by an atom is assumed to be elastic (no loss in momentum), the magnitudes of the wave vectors are the same. Thus, s – s0 becomes the base of an isosceles triangle with equal sides of Is0I and IsI. We have:

\left | \textbf{\textit{s}}-\textbf{\textit{s}}_0\right |=\left | \textbf{\textit{s}}\right |sin\theta +\left | \textbf{\textit{s}}_0\right |sin\theta

Since \left | \textbf{\textit{s}}\right |=\left | \textbf{\textit{s}}_0\right |=\frac{1}{\lambda } (see previous article)

\left | \textbf{\textit{s}}-\textbf{\textit{s}}_0\right |=\frac{2sin\theta }{\lambda }\; \; \; \; \; \; \; (23)

To find an expression for dnh,nk,nl (see above diagram), we determine the scalar projection of a/h (the vector linking two plane-intercept points on the a-axis) on the vector s – s0 to give:

d_{nh,nk,nl}=\frac{\textbf{\textit{a}}}{h}\cdot \frac{\textbf{\textit{s}}-\textbf{\textit{s}}_0}{\left | \textbf{\textit{s}}-\textbf{\textit{s}}_0\right |}

Dividing both sides of eq20 by h, we have \frac{\textbf{\textit{a}}}{h}\cdot ( \textbf{\textit{s}}-\textbf{\textit{s}}_0) =1. So,

d_{nh,nk,nl}=\frac{1}{\left | \textbf{\textit{s}}-\textbf{\textit{s}}_0\right |}\; \; \; \; \; \; \; (24)

Substitute eq23 in eq24, we have the Bragg equation:

2d_{nh,nk,nl}\: sin\theta =\lambda

Since the Laue equations are equal to the Bragg equation if we denote the angle between the wave vectors s and s0  by 2θ, the Bragg equation is a specific form of the Laue equations with the condition of 2θ imposed.

 

next article: reciprocal lattice
Previous article: Laue equations
Content page of X-ray crystallography
Content page of advanced chemistry
Main content page