Category: Advanced Chemistry

Posted on

Constitutive relation: hypothesis 4

The fourth constitutive relation hypothesis states that the properties of a fluid are isotropic, i.e., independent of direction.

For example, an object moving in a fluid encounters the same resistance regardless of the direction of movement. From eq11, the properties of a fluid are described by βijkl , which must be isotropic.

As described in the articles on tensors, the general form of a fourth-order isotropic tensor is:

\beta_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu\delta_{ik}\delta_{jl}+\nu\delta_{il}\delta_{jk}\; \; \; \; \; \; \; (12)

From eq11 of the previous article, τij is symmetric. Therefore, τij =τji and

\beta_{ijkl}=\beta_{jikl}

Substitute eq12 in the above equation, we have

\left ( \mu-\nu \right )\delta_{ik}\delta_{jl}=\left ( \mu-\nu \right )\delta_{il}\delta_{jk}

\mu=\nu\; \; \; \; \; \; \; (13)

Substitute eq13 in eq12,

\beta_{ijkl}=\lambda\delta_{ij}\delta_{kl}+\mu\left (\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right )\; \; \; \; \; \; \; (14)

Substitute eq14 in eq11,

\tau_{ij}=\frac{1}{2}\lambda\delta_{ij}\delta_{kl}\left ( \frac{\partial u_k}{\partial x_l}+\frac{\partial u_l}{\partial x_k}\right ) +\frac{1}{2}\mu\left (\delta_{ik}\delta_{jl}+\delta_{il}\delta_{jk}\right )\left ( \frac{\partial u_k}{\partial x_l}+\frac{\partial u_l}{\partial x_k}\right )

\tau_{ij}=\frac{1}{2}\lambda\delta_{ij}\left ( \frac{\partial u_k}{\partial x_k}+\frac{\partial u_k}{\partial x_k}\right ) +\frac{1}{2}\mu\left (\delta_{ik}\frac{\partial u_k}{\partial x_j}+\delta_{ik}\frac{\partial u_j}{\partial x_k} +\delta_{il}\frac{\partial u_j}{\partial x_l}+\delta_{il}\frac{\partial u_l}{\partial x_j} \right )

\tau_{ij}=\lambda\delta_{ij} \frac{\partial u_k}{\partial x_k}+\mu\left (\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right )\; \; \; \; \; \; \; (15)

Substitute eq15 in eq3,

\sigma_{ij}=-p\delta_{ij}+\lambda\delta_{ij} \frac{\partial u_k}{\partial x_k}+\mu\left (\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right )\; \; \; \; \; \; \; (16)

The values of μ and λ can only be determined through experiments. μ is known as the shear viscosity of the fluid while λ is the volume viscosity of the fluid, which is zero for an incompressible fluid. Eq16 becomes:

\sigma_{ij}=-p\delta_{ij}+\mu\left (\frac{\partial u_i}{\partial x_j}+\frac{\partial u_j}{\partial x_i} \right )\; \; \; \; \; \; \; (17)

 

Question

Show that μ in eq16 and eq7 are the same.

Answer

Consider the flow of an incompressible fluid where the flow velocity components are u = u(y), v = 0 and w = 0. Using eq16, the stress components are:

\sigma_{11}=\sigma_{22}=\sigma_{33}=-p

\sigma_{12}=\sigma_{21}=\mu\frac{du}{dy}

with the remaining components equal to zero.

Clearly, the components of stress in this case include the pressure p and the shear stress  \mu\frac{du}{dy}, which is the result of the derivation of Newton’s law of viscosity, eq7. Hence, μ in eq16 and eq7 are the same.

 

 

Previous article: hypothesis 3
Content page of constitutive relation
Content page of advanced chemistry
Main content page
Posted on

Constitutive relation: hypothesis 2

The second hypothesis, formulating the constitutive relation for an incompressible and viscous Newtonian fluid, states that the shear stress tensor is a linear function of the deformation tensor.

Consider a shear stress τ that is coplanar to a cross section of a Newtonian fluid ABCD as shown in the diagram below.

The fluid is continuously deformed by the force with different layers of the fluid having different velocities, for example, u(y) at point D and u(y + Δy) at point A. At time t, the body of fluid moves from ABCD to A’B’CD. Let’s define the shear strain in the fluid as \gamma=\frac{\Delta x}{\Delta y} where

\frac{\Delta x}{\Delta y}=\frac{\left [ u\left ( y+\Delta y \right )-u\left ( y \right ) \right ]\Delta t}{\Delta y}

If Δy → 0,

\gamma=\frac{du}{dy}\Delta t

If t is small,

\frac{d\gamma}{dt}=\frac{du}{dy}\; \; \; \; \; \; \; (5)

We call \frac{d\gamma}{dt} the shear strain rate and \frac{du}{dy} the velocity gradient of the fluid. Clearly, the greater the shear stress applied, the higher the value of the shear strain rate, i.e.:

\tau \propto \frac{d\gamma}{dt} = \mu\frac{d\gamma}{dt}\; \; \; \; \; \; \; (6)

where μ is the constant of proportionality, which is attributed to the resistance of the fluid to the shear stress.

We call this resistance, the viscosity of the fluid. Substituting eq5 in eq6, we have the Newton’s law of viscosity,

\tau=\mu\frac{du}{dy}\; \; \; \; \; \; \; (7)

When extended to three dimensions (see above diagram), eq7 becomes:

\tau_{yx}=\mu\frac{du_x}{dy}\; \; \; \; \; \; \; (7a)

From eq4 of the previous article, shear stress is a second-order tensor with nine components τij . Hence, the general form of eq7 is:

\tau_{ij}=\mu\frac{\partial u_j}{\partial x_i}

where for a Cartesian system, i = 1, 2, 3, j = 1, 2, 3 and x1  = x, x2  = y, x3  = z.

\frac{\partial u_j}{\partial x_i} can therefore be expressed as a second-order tensor (deformation tensor) of nine components. Since each of the nine components of τij is linearly proportional to each of the nine components of \frac{\partial u_j}{\partial x_i} , there are a total of 92 = 81 proportionality constants. In other words,

\tau_{ij}=\alpha_{ijkl}\frac{\partial u_k}{\partial x_l}\; \; \; \; \; \; \; (8)

where \alpha_{ijkl} is a fourth-order tensor of eighty one components, which can be representation by a 3x3x3x3 matrix for a Cartesian system. This satisfies hypothesis 2.

The fourth-order tensor, which is composed of the proportionality constants, characterises the properties of the fluid.

 

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

Isotropic tensors

An isotropic Cartesian tensor is one where its components are identical in any orthogonal Cartesian system.

An isotropic property is one that is independent of direction, e.g. thermal expansion of a solid, and therefore has a quantity that is independent of the reference frame. An isotropic property can therefore be described by a Cartesian tensor like T'_{j_1j_2...j_n} in eq15. 

To evaluate the properties of an isotropic tensor, we begin by considering a generic fourth-order tensor T_{jlnp}, whose transformation is expressed by eq14:

T'_{ikmo}=a_{ij}a_{kl}a_{mn}a_{op}T_{jlnp}\; \; \; \; \; \; \; (16)

Since j, l, n and p, each of which ranges from 1 to 3, the 81 components of T'_{ikmo} are classified in four groups as follows:

    1. Components with four equal subscripts, i.e. T'_{1111}, T'_{2222} and T'_{3333}.
    2. Components with three equal subscripts, e.g. T'_{1211}, T'_{2223}, etc.
    3. Components with two different pairs of equal subscripts, e.g. T'_{1212}, T'_{3322}, T'_{1331}, etc.
    4. Components with one pair of equal subscripts, e.g. T'_{1213}, T'_{2213}, etc.

Next, consider the transformation of T_{jlnp} to T'_{ikmo} as a rotation about the x3-axis of \theta=180^o where x3 coincides with x’3 (x’3 is equal to x3 and is perpendicular to the plane of the page). From eq12, the transformation matrix is:

a_{ij}\: or\: a_{kl}\: or\: a_{mn}\: or\: a_{op}=\begin{pmatrix} cos\theta &sin\theta &0 \\ -sin\theta &cos\theta &0 \\ 0 &0 &1 \end{pmatrix}=\begin{pmatrix} -1 &0 &0 \\ 0 &-1 &0 \\ 0 &0 &1 \end{pmatrix}

Therefore, a_{11}=a_{22}=-1 and a_{33}=1, with the rest of the components of the transformation matrix equal to zero. From eq16, a component of group 2 of the above classification is:

T'_{1113}=a_{1j}a_{1l}a_{1n}a_{3p}T_{jlnp}=a_{11}a_{11}a_{11}a_{31}T_{1111}+a_{11}a_{11}a_{11}a_{32}T_{1112}

+a_{11}a_{11}a_{11}a_{33}T_{1113}+...\; \; \; \; \; \; \; (17)

Substituting a_{11}=a_{22}=-1, a_{33}=1 and a_{i\neq j}=0 in eq17, we have:

T'_{1113}=-T_{1113}\; \; \; \; \; \; \; (18)

If T_{jlnp} is an isotropic tensor, T'_{1113}=T_{1113} (since an isotropic Cartesian tensor is one where its components are identical in any orthogonal Cartesian system). Hence, the only way to satisfy eq18 is for T'_{1113}=T_{1113}=0. Otherwise, the value of T'_{1113} varies with the value of T_{1113}, e.g. when T_{1113}=-7, T'_{1113}=7.

Another component of group 2 of the above classification is T'_{1211} where after expanding and substituting with a_{11}=a_{22}=-1, a_{33}=1 and a_{i\neq j}=0,

T'_{1211}=T_{1211}\; \; \; \; \; \; \; (20)

Eq20 shows that T'_{1211}=T_{1211} but does not indicate the value of either component. If the rotation is instead made about the x2-axis at \theta=180^o, the transformation matrix is:

\begin{pmatrix} cos\theta & 0 &sin\theta \\ 0 & 1 &0 \\ -sin\theta &0 &cos\theta \end{pmatrix}=\begin{pmatrix} -1 &0 &0 \\ 0 &1 &0 \\ 0 &0 &-1 \end{pmatrix}

We have a_{11}=a_{33}=-1,a_{22}=1, and a_{i\neq j}=0, giving T'_{1211}=-T_{1211}. If T_{jlnp} is an isotropic tensor,

T'_{1211}=T_{1211}=0\; \; \; \; \; \; \; (22)

Repeating the above steps and assuming that T_{jlnp} is an isotropic tensor, we find that all the components of T'_{ikmo} in group 2 and group 4 are zero, which means that all the components of T_{jlnp} in group 2 and group 4 are zero.

We shall now analyse the components of group 1 by a rotation about the x3-axis at \theta=90^o with the transformation matrix:

\begin{pmatrix} cos\theta &sin\theta &0 \\ -sin\theta &cos\theta &0 \\ 0 &0 &1 \end{pmatrix}=\begin{pmatrix} 0 &1 &0 \\ -1 &0 &0 \\ 0 &0 &1 \end{pmatrix}

We have T'_{1111}=T_{2222},T'_{2222}=T_{1111},T'_{3333}=T_{3333}.

A rotation about the x2-axis at \theta=90^o involves the transformation matrix:

\begin{pmatrix} cos\theta & 0&sin\theta \\ 0 &1 &0 \\ -sin\theta &0 & cos\theta\end{pmatrix}=\begin{pmatrix} 0 &0 &1 \\ 0 &1 &0 \\ -1 &0 &0 \end{pmatrix}

and gives T'_{1111}=T_{3333},T'_{2222}=T_{2222},T'_{3333}=T_{1111}.

Comparing the six equations T'_{1111}=T_{2222},T'_{2222}=T_{1111},T'_{3333}=T_{3333}  and T'_{1111}=T_{3333},T'_{2222}=T_{2222},T'_{3333}=T_{1111}, all the components of T'_{ikmo} in group 1 are equal to one another, which means that all the components of T_{jlnp} in group 1 are equal to one another.

Finally, let’s look at the components in group 3. Under a rotation about the x3-axis at \theta=90^o, we have T'_{1122}=T_{2211},T'_{2233}=T_{1133},T'_{3311}=T_{3322}.

A rotation about the x2-axis at \theta=90^o gives T'_{1122}=T_{3322},T'_{2233}=T_{2211},T'_{3311}=T_{1133}.

Comparing the six equations T'_{1122}=T_{2211},T'_{2233}=T_{1133},T'_{3311}=T_{3322}  and T'_{1122}=T_{3322},T'_{2233}=T_{2211},T'_{3311}=T_{1133} , components in group 3 of T'_{ikmo} (which is equal to the components in group 3 of T_{jlnp}) with i = k and m = o are equal to one another. Repeating the above steps, we can show that components in group 3 with i = m and k = o are equal to one another and those with i = o and k = m are equal to one another.

Therefore, if T_{jlnp} is an isotropic tensor:

i) components of T'_{ikmo} with i = k = m = o are equal and therefore T_{jlnp} with j = l = n = p are equal

ii) components of T'_{ikmo} with i = k ≠ m = o are equal and therefore T_{jlnp} with j = l ≠ n = p are equal

iii) components of T'_{ikmo} with i = m ≠ k = o are equal and therefore T_{jlnp} with j = n ≠ l = p are equal

iv) components of T'_{ikmo} with i = o ≠ k = m are equal and therefore T_{jlnp} with j = p ≠ l = n are equal

v) all other components not mentioned are equal to zero

Note that components of parts ii, iii and iv belong to group3 of the isotropic tensor. Now, consider a case of T_{jlnp} where components of T_{jlnp} in iii and iv are set to zero (T_{jlnp} is still an isotropic tensor but with more components equal to zero). A rotation about the x3-axis at \theta=45^o gives the transformation matrix:

\begin{pmatrix} cos\theta &sin\theta &0 \\ -sin\theta &cos\theta &0 \\ 0 & 0 &1 \end{pmatrix}=\begin{pmatrix} \frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} & 0\\ -\frac{\sqrt{2}}{2} & \frac{\sqrt{2}}{2} &0 \\ 0 &0 &1 \end{pmatrix}

Since components in part iii and part iv are set to zero and the conditions of part v applies, we have T'_{1111}=\frac{1}{4}T_{1111}+\frac{1}{4}T_{1122}+\frac{1}{4}T_{2211}+\frac{1}{4}T_{2222} . Furthermore, part i states that T_{1111}=T_{2222} and part ii states that T_{1122}=T_{2211} . So,

T'_{1111}=\frac{1}{2}\left ( T_{1111}+T_{1122}\right )\; \; \; \; \; \; \; (23)

If T_{jlnp} is an isotropic tensor,  T'_{1111}=T_{1111} and therefore, T_{1111}=T_{1122} . This means that all non-zero components of this particular T_{jlnp} have the same value, say, λ, where:

T_{jlnp}=\lambda\delta_{jl}\delta_{np}\; \; \; \; \; \; \; (24)

where \delta_{jl} and \delta_{np} are Kronecker deltas.

Using the same logic by considering the cases of T_{jlnp} where components of T_{jlnp} in part ii and part iv are set to zero and where components of T_{jlnp} in part ii and part iii are set to zero, we have:

T_{jlnp}=\mu\delta_{jn}\delta_{lp}\; \; \; \; \; \; \; (25)

and

T_{jlnp}=\nu\delta_{jp}\delta_{ln}\; \; \; \; \; \; \; (26)

respectively.

Therefore, we can write the general form of the fourth-order isotropic tensor as:

T_{jlnp}=\lambda\delta_{jl}\delta_{np}+\mu\delta_{jn}\delta_{lp}+\nu\delta_{jp}\delta_{ln}\; \; \; \; \; \; \; (27)

 

Question

Show that T_{jlnp} in eq27 continues to satisfy parts i to v, all of which are needed to define a fourth-order isotropic tensor.

Answer

In eq27,

i) the components of T_{jlnp} with j = l = n = p have values of \lambda+\mu+\nu and are therefore still equal to one another;

ii) the components of T_{jlnp} with j = l ≠ n = p have values of \lambda and are therefore still equal to one another;

iii) the components of T_{jlnp} with j = n ≠ l = p have values of \mu and are therefore still equal to one another;

iv) the components of T_{jlnp} with j = p ≠ l = n have values of \nu and are therefore still equal to one another;

v) all other components not mentioned are equal to zero.

Therefore, the general form of the fourth-order isotropic tensor is valid.

 

 

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

Tensor and tensor transformation

A tensor is an array of numbers or functions that can be used to describe properties of a body, such as the scalar for temperature, the vector \begin{pmatrix} a_1\\a_2 \\a_3 \end{pmatrix}  for velocity, and the matrix \begin{pmatrix} \sigma_{11} &\sigma_{12}&\sigma_{13}\\ \sigma_{21} &\sigma_{22}&\sigma_{23}\\ \sigma_{31} &\sigma_{32}& \sigma_{33}\end{pmatrix}  for stress.

Tensors are categorised by their ranks (also known as order), where

    1. A zeroth-order tensor is an array representing a quantity that only has magnitude and no direction, i.e. a scalar.
    2. A first-order tensor is an array representing a quantity that has magnitude and a direction, i.e. a vector.
    3. A second-order tensor is an array representing a quantity that has magnitude and two directions.

The above implies that the number of elements of a tensor is 3n, where n is the order of the tensor.

Next, we shall investigate how the elements of a tensor, particularly a 2nd-order tensor, transform between two sets of coordinates. From an earlier article, uj, which is a set of three quantities, is known as a first order tensor (i.e. a vector). Consider the product of two first order tensors uj and vl , whose elements with respect to another set of axes are given by:

u'_iv'_k=\left ( a_{ij}u_j \right )\left ( a_{kl}v_l \right )\; \; \; \; \; \; \; \; 12a

Since each term on the RHS of eq12a is a scalar, we can rewrite eq12a as

u'_iv'_k=a_{ij}a_{kl}u_jv_l=\sum_{j}\sum_{i}a_{ij}a_{kl}u_jv_l\; \; \; \; \; \; \; \; 12b

There are 9 possible products of u'_iv'_k because i = 1,2,3 and k = 1,2,3. In other words, u'_iv'_k is a set of 9 elements, each being a sum of nine quantities , e.g.

u'_1v'_1=a_{11}a_{11}u_1v_1+a_{11}a_{12}u_1v_2+a_{11}a_{13}u_1v_3

+a_{12}a_{11}u_2v_1+a_{12}a_{12}u_2v_2+a_{12}a_{13}u_2v_3

+a_{13}a_{11}u_3v_1+a_{13}a_{12}u_3v_2+a_{13}a_{13}u_3v_3

We can therefore represent the elements of u'_iv'_k by a 3×3 matrix \begin{pmatrix} u'_1v'_1 &u'_1v'_2&u'_1v'_3\\ u'_2v'_1&u'_2v'_2&u'_2v'_3\\ u'_3v'_1&u'_3v'_2&u'_3v'_3 \end{pmatrix}. Similarly, u_jv_l in eq12b is a multiplication of two first order tensors and can be represented by a 3×3 matrix \begin{pmatrix} u_1v_1 &u_1v_2&u_1v_3\\ u_2v_1&u_2v_2&u_2v_3\\ u_3v_1&u_3v_2&u_3v_3 \end{pmatrix}. This results in:

T'_{ik}=a_{ij}a_{kl}T_{jl}\; \; \; \; \; \; \; (13)

Eq13 indicates that each element of u'_1v'_1 in one reference frame is a result of the transformation of u_jv_l in a different reference frame by a_{ij}a_{kl} . In general, a set of nine quantities T_{jl} (as is the case of a second-order tensor) with reference to an orthogonal set of axes is transformed to another set of nine quantities T'_{ik} with reference to another orthogonal set of axes by the transformation matrix a_{ij}a_{kl}.

Using the same logic, a third order tensor transforms as follows:

u'_iv'_kw'_m=\left ( a_{ij}u_j \right )\left ( a_{kl}u_l \right )\left ( a_{mn}u_n \right )=a_{ij}a_{kl}a_{mn}u_jv_lw_n

T'_{ikm}=a_{ij}a_{kl}a_{mn}T_{jln}

and a fourth order tensor:

u'_iv'_kw'_mx'_o=\left ( a_{ij}u_j \right )\left ( a_{kl}u_l \right )\left ( a_{mn}u_n \right )\left ( a_{op}x_p \right )=a_{ij}a_{kl}a_{mn}a_{op}u_jv_lw_nx_p

T'_{ikmo}=a_{ij}a_{kl}a_{mn}a_{op}T_{jlnp}\; \; \; \; \; \; \; (14)

In general, a tensor of rank n is a quantity that transforms from one set of orthogonal axes to another set of orthogonal axes according to:

T'_{j_1,j_2...j_n}=a_{j_1k_1}a_{j_2k_2}...a_{j_nk_n}T_{k_1k_2...k_n}\; \; \; \; \; \; \; (15)

 

previous article: orthogonal transformation
Next article: isotropic tensors
Content page of tensors
Content page of advanced chemistry
Main content page
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