Tensor Operations in OpenFOAM


Tensor operations are frequently used in turbulence modeling because the strain rate and vorticity tensors can be expressed in terms of the velocity gradient tensor.

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

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

OpenFOAM Version
OpenFOAM-dev

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

Author: fumiya

CFD engineer in Japan

Leave a Reply

Your email address will not be published. Required fields are marked *