$Q$ 准则是一种常用的漩涡识别方法，由 Hunt 等1提出。要了解 $Q$ 的定义，首先要知道什么是张量的不变量。

张量的不变量

对于 $3\times3$ 的张量 $\mathbf A$，其特征方程可以表示为如下形式

$$ \begin{equation} \lambda^3 - I_1 \lambda^2 + I_2 \lambda - I_3 = 0 \end{equation} $$

其中，$\lambda$ 是 $\mathbf A$ 的特征值，$I_1$、$I_2$、$I_3$ 分别称为 $\mathbf A$ 的第一、第二、第三不变量，定义如下

$$ \begin{equation} I_1 = \mathrm{tr} (\mathbf A) \end{equation} $$

$$ \begin{equation} I_2 = \frac{1}{2} \left[ ( \mathrm{tr}(\mathbf A))^2 - \mathrm{tr}(\mathbf A^2 ) \right] \label{eq:2nd_invariant} \end{equation} $$

$$ \begin{equation} I_3 = \mathrm{det}(\mathbf A) \end{equation} $$

其中， $\mathrm{tr} (\mathbf A)$ 表示 $\mathbf A$ 的迹（trace），$\mathrm{det} (\mathbf A)$ 表示 $\mathbf A$ 的行列式（determinant）。

张量的迹定义如下：

$$ \begin{equation} \mathrm{tr} (\mathbf A) = A_{11} + A_{22} + A_{33} \end{equation} $$

若张量为矢量的梯度（如 $\nabla \mathbf u$），则有如下关系

$$ \begin{equation} \mathrm{tr} (\nabla \mathbf u) = \frac{\partial u}{\partial x} + \frac{\partial v}{\partial y} + \frac{\partial w}{\partial z} = \nabla \cdot \mathbf u \end{equation} $$

Q 的定义

$Q$ 的定义为速度梯度的第二不变量（second invariant of velocity gradient）。将方程 $\eqref{eq:2nd_invariant}$ 中的 $\mathbf A$ 替换为速度梯度 $\nabla \mathbf u$，就得到了 $Q$

$$ \begin{equation} Q = \frac{1}{2} \left[ ( \mathrm{tr}(\nabla \mathbf u))^2 - \mathrm{tr}((\nabla \mathbf u)^2 ) \right] \label{eq:q_criterion} \end{equation} $$

式 $\eqref{eq:q_criterion}$ 为 OpenFOAM 中计算 $Q$ 时所用的公式。

对于不可压缩流体，上式可以进一步化简。根据连续性方程，有

$$ \begin{equation} \mathrm{tr} (\nabla \mathbf u) = \nabla \cdot \mathbf u = 0 \end{equation} $$

式 $\eqref{eq:q_criterion}$ 等号右边第一项为 0，变成了

$$ \begin{equation} Q = - \frac{1}{2} \mathrm{tr}((\nabla \mathbf u)^2 ) \end{equation} $$

将速度梯度分解成两部分

$$ \begin{equation} \nabla \mathbf u = \underbrace{\frac{1}{2}\left(\nabla \mathbf u + (\nabla \mathbf u)^T \right)}_ % workaround for hugo {\text{对称张量}}+\underbrace{\frac{1}{2}\left(\nabla \mathbf u - (\nabla \mathbf u)^T \right)}_{\text{反对称张量}} \end{equation} $$

上式中的对称张量记为 $\mathbf S$，被称为应变张量 （strain rate tensor）；反对称张量记为 $\mathbf \Omega$，被称为旋转张量（spin tensor）

$$ \begin{equation} \mathbf S =\frac{1}{2}\left(\nabla \mathbf u + (\nabla \mathbf u)^T \right) \end{equation} $$

$$ \begin{equation} \mathbf \Omega =\frac{1}{2}\left(\nabla \mathbf u - (\nabla \mathbf u)^T \right) \end{equation} $$

$Q$ 可写作

$$ \begin{equation} \begin{aligned} Q = -\frac{1}{2}\mathrm{tr}((\nabla \mathbf u)^2 )\mathbf \newline & = -\frac{1}{2}\mathrm{tr}(\mathbf S^2 + 2 \mathbf S \mathbf \Omega + \mathbf \Omega^2) \newline & = -\frac{1}{2}\left[ \mathrm{tr}(\mathbf S^2) + \mathrm{tr} (\mathbf S \mathbf \Omega) +\mathrm{tr}(\mathbf \Omega^2) \right] \newline \end{aligned} \end{equation} $$

由于 $\mathbf S$ 和 $\mathbf \Omega$ 分别是对称张量和反对称张量，因此具有以下性质

$$ \begin{equation} \mathbf S = \mathbf S^T \end{equation} $$

$$ \begin{equation} \mathbf \Omega = - \mathbf \Omega^T \end{equation} $$

所以有

$$ \begin{equation} \mathrm{tr}(\mathbf S^2) = \mathrm{tr}(\mathbf S^T \mathbf S) \end{equation} $$

$$ \begin{equation} \mathrm{tr}(\mathbf \Omega^2) = - \mathrm{tr}(\mathbf \Omega^T \mathbf \Omega) \end{equation} $$

且对于 $\mathbf S$ 和 $\mathbf \Omega$，有如下公式

$$ \begin{equation} \begin{aligned} \mathrm{tr} (\mathbf S \mathbf \Omega) & = \mathrm{tr} \left( \frac{1}{2}\left( \nabla\mathbf u + \left( \nabla \mathbf u \right)^T \right) \frac{1}{2}\left( \nabla\mathbf u - \left( \nabla \mathbf u \right)^T \right) \right) \newline & = \mathrm{tr} \left( \frac{1}{4} \left( (\nabla \mathbf u)^2 - ((\nabla \mathbf u)^T)^2 \right) \right) \newline & = 0 \end{aligned} \end{equation} $$

于是 $Q$ 最终可以写作

$$ \begin{equation} Q = -\frac{1}{2} \left[\mathrm{tr}(\mathbf S^T \mathbf S) - \mathrm{tr}(\mathbf \Omega^T \mathbf \Omega) \right] = \frac{1}{2}\left[ || \mathbf \Omega ||_F^2 - || \mathbf S ||_F^2 \right] \label{eq:ico_q_criterion} \end{equation} $$

其中 $||\quad||_F$ 为 F-范数（Frobenius norm），定义为

$$ \begin{equation} ||\mathbf A||_F = \sqrt{\mathrm{tr}(\mathbf A^T \mathbf A)} \end{equation} $$

通常将 $\eqref{eq:ico_q_criterion}$ 简写为

$$ \begin{equation} Q=\frac{1}{2}\left[ || \mathbf \Omega ||^2 - || \mathbf S ||^2 \right] \label{eq:ico_q_criterion_simp} \end{equation} $$

式 $\eqref{eq:ico_q_criterion_simp}$ 是文献中通常采用的形式。$Q>0$ 表示旋转程度大于变形程度，认为该区域存在漩涡。

代码实现

OpenFOAM 3.0.x

OpenFOAM 3.0.x 计算 $Q$ 的代码有两份：一份是 utility ，编译后生成可执行文件，可以直接用Q命令调用；另一份是 function object，编译后生成动态链接库文件，通过 function objects 机制调用。二者的代码分别位于：

  • applications/utilities/postProcessing/velocityField/Q/
  • src/postProcessing/functionObjects/utilities/Q/

utility 中的源码中提供了两种计算方式，第一种采用公式 $\eqref{eq:q_criterion}$ 计算，代码如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14

        volTensorField gradU(fvc::grad(U));

        volScalarField Q
        (
            IOobject
            (
                "Q",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            0.5*(sqr(tr(gradU)) - tr(((gradU)&(gradU))))
        );

第二种采用公式 $\eqref{eq:ico_q_criterion_simp}$ 计算，是不可压缩流的计算方式，代码如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

        // This is a second way of calculating Q, that delivers results
        // very close, but not identical to the first approach.

        volSymmTensorField S(symm(gradU));  // symmetric part of tensor
        volTensorField W(skew(gradU));      // anti-symmetric part

        volScalarField SS(S && S);
        volScalarField WW(W && W);
        volScalarField Q
        (
            IOobject
            (
                "Q",
                runTime.timeName(),
                mesh,
                IOobject::NO_READ,
                IOobject::NO_WRITE
            ),
            0.5*(WW - SS)
        );

第二种方式在代码中是被注释掉的，默认采用第一种方式计算。

function object 中的源码采用公式 $\eqref{eq:q_criterion}$ 计算，代码如下：

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19
20

void Foam::Q::execute()
{
    if (active_)
    {
        const fvMesh& mesh = refCast<const fvMesh>(obr_);

        const volVectorField& U =
            mesh.lookupObject<volVectorField>(UName_);

        const volTensorField gradU(fvc::grad(U));

        volScalarField& Q =
            const_cast<volScalarField&>
            (
                mesh.lookupObject<volScalarField>(type())
            );

        Q = 0.5*(sqr(tr(gradU)) - tr(((gradU) & (gradU))));
    }
}

OpenFOAM 4.x 及后续版本

OpenFOAM 4.x 将大部分 utilities 合并到 functionObjects ，并将 functionObjects 框架重写，计算 $Q$ 的代码位于：

  • src/functionObjects/field/Q/

采用公式 $\eqref{eq:q_criterion}$ 计算，代码如下

 1
 2
 3
 4
 5
 6
 7
 8
 9
10
11
12
13
14
15
16
17
18
19

bool Foam::functionObjects::Q::calc()
{
    if (foundObject<volVectorField>(fieldName_))
    {
        const volVectorField& U = lookupObject<volVectorField>(fieldName_);
        const tmp<volTensorField> tgradU(fvc::grad(U));
        const volTensorField& gradU = tgradU();

        return store
        (
            resultName_,
            0.5*(sqr(tr(gradU)) - tr(((gradU) & (gradU))))
        );
    }
    else
    {
        return false;
    }
}

参考资料

  1. Hunt, J. C. R., Wray, A. and Moin, P. (1988). Eddies, Stream, and Convergence Zones in Turbulent Flows. Center for Turbulence Research. Report CTR-S88. ↩︎