Tensor Operations in OpenFOAM

Last Updated: May 5, 2019

The gradient of the velocity field is a strain-rate tensor field, that is, a second rank tensor field. It appears in the diffusion term of the Navier-Stokes equation.

A second rank tensor has nine components and can be expressed as a 3×3 matrix as shown in the above image.

In this blog post, I will pick out some typical tensor operations and give brief explanations of them with some usage examples in OpenFOAM.

Keywords
strain rate tensor, vorticity tensor, Q-criterion, Hodge dual

Gradient of a Vector Field | fvc::grad(u)

The gradient of a velocity vector \(\boldsymbol{u}\) returns a velocity gradient tensor (second rank tensor).
\begin{eqnarray}
\nabla \boldsymbol{u} &\equiv& \partial_{i} u_j \\
&=& \left(
\begin{matrix}
\partial u_1/\partial x_1 & \partial u_2/\partial x_1 & \partial u_3/\partial x_1 \\
\partial u_1/\partial x_2 & \partial u_2/\partial x_2 & \partial u_3/\partial x_2 \\
\partial u_1/\partial x_3 & \partial u_2/\partial x_3 & \partial u_3/\partial x_3
\end{matrix} \right)
\end{eqnarray}

Symmetric Part of a Second Rank Tensor | symm(T)

\begin{eqnarray}
{\rm symm}(\boldsymbol{T}) &\equiv& \frac{1}{2} (\boldsymbol{T} + \boldsymbol{T}^T) \\
&=& \frac{1}{2} \left(
\begin{matrix}
2T_{11} & T_{12} + T_{21} & T_{13} + T_{31} \\
T_{21} + T_{12} & 2T_{22} & T_{23} + T_{32} \\
T_{31} + T_{13} & T_{32} + T_{23} & 2T_{33}
\end{matrix} \right)
\end{eqnarray}

Twice the Symmetric Part of a Second Rank Tensor | twoSymm(T)

\begin{eqnarray}
{\rm twoSymm}(\boldsymbol{T}) &\equiv& \boldsymbol{T} + \boldsymbol{T}^T \\
&=& \left(
\begin{matrix}
2T_{11} & T_{12} + T_{21} & T_{13} + T_{31} \\
T_{21} + T_{12} & 2T_{22} & T_{23} + T_{32} \\
T_{31} + T_{13} & T_{32} + T_{23} & 2T_{33}
\end{matrix} \right)
\end{eqnarray}

Skew-symmetric Part of a Second Rank Tensor | skew(T)

\begin{eqnarray}
{\rm skew}(\boldsymbol{T}) &\equiv& \frac{1}{2} (\boldsymbol{T} – \boldsymbol{T}^T) \\
&=& \frac{1}{2} \left(
\begin{matrix}
0 & T_{12} – T_{21} & T_{13} – T_{31} \\
T_{21} – T_{12} & 0 & T_{23} – T_{32} \\
T_{31} – T_{13} & T_{32} – T_{23} & 0
\end{matrix} \right)
\end{eqnarray}

The symm and skew operations of the velocity gradient tensor field frequently appear in the source code. The following is a typical example.

In the above code, the symmetric and antisymmetric parts of the velocity gradient tensor \(\partial u_j/\partial x_i\) are defined as follows:
\begin{eqnarray}
S_{ij} &=& \frac{1}{2} \left( \frac{\partial u_j}{\partial x_i} + \frac{\partial u_i}{\partial x_j} \right), \\
\Omega_{ij} &=& \frac{1}{2} \left( \frac{\partial u_j}{\partial x_i} – \frac{\partial u_i}{\partial x_j} \right),
\end{eqnarray}
where \(S_{ij}\) is the strain rate tensor and \(-\Omega_{ij}\) is the vorticity (spin) tensor.

Hodge Dual | *T

\begin{equation}
*T = \left( T_{23},\;-T_{13},\;T_{12} \right)
\end{equation}

The vorticity vector \(\boldsymbol{\omega}\) is calculated as the Hodge dual of the skew-symmetric part of the velocity gradient tensor.
\begin{eqnarray}
\boldsymbol{\omega} &=& 2 \times \left( * \Omega_{ij} \right) \\
&=& 2 \times \left( \Omega_{23},\;-\Omega_{13},\;\Omega_{12} \right) \\
&=& \left( \frac{\partial u_3}{\partial x_2}-\frac{\partial u_2}{\partial x_3},\;\frac{\partial u_1}{\partial x_3}-\frac{\partial u_3}{\partial x_1},\;\frac{\partial u_2}{\partial x_1}-\frac{\partial u_1}{\partial x_2} \right)
\end{eqnarray}

The operator * in front of the skew is the Hodge dual operator.

Inner Product of Two Second Rank Tensors | T & S

\begin{equation}
P_{ij} = T_{ik} S_{kj}
\end{equation}

Double Inner Product of Two Second Rank Tensors | T && S

\begin{align}
s = T_{ij}S_{ij} &= T_{11}S_{11} + T_{12}S_{12} + T_{13}S_{13} \\
&+ T_{21}S_{21} + T_{22}S_{22} + T_{23}S_{23} \\
&+ T_{31}S_{31} + T_{32}S_{32} + T_{33}S_{33}
\end{align}

Trace of a Second Rank Tensor | tr(T)

\begin{equation}
{\rm tr}(\boldsymbol{T}) \equiv T_{11} + T_{22} + T_{33} = \boldsymbol{T} {\rm \&}{\rm \&} \boldsymbol{I}
\end{equation}

The Q-criterion can be used to identify vortex cores:
\begin{align}
Q &= \frac{1}{2} \left[ (\nabla \cdot \boldsymbol{u})^2 – {\rm tr}(\nabla \boldsymbol{u}^2)\right] \\
&= \frac{1}{2} \left[ (\nabla \cdot \boldsymbol{u})^2 + \Omega_{ij}\Omega_{ij} – S_{ij}S_{ij} \right].
\end{align}
For incompressible flows, it can be simplified as follows:
\begin{equation}
Q = \frac{1}{2} \left[ \Omega_{ij}\Omega_{ij} – S_{ij}S_{ij} \right].
\end{equation}

References

[1] OpenFOAM Programmer’s Guide
[2] CFD Direct | Tensor Mathematics