在 OpenFOAM4.1 中显示三维涡结构,核心方法是结合使用软件自带的 Q-criterion(Q准则)或 Lambda2(λ₂准则)等旋涡识别方法,并在后处理软件 ParaView 中进行可视化。

以下是具体的操作步骤与实现方法:

1. 在 OpenFOAM 中提取涡旋数据

您可以在计算过程中或计算后通过添加函数对象(Function Objects)来提取涡量、Q准则或 Lambda2。
https://www.topcfd.cn/12128/

  • 打开算例目录下的 system/controlDict 文件。
  • 在文件末尾添加以下 functions 字典:
1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
functions
{
Q1
{
type Q;
libs ("libfieldFunctionObjects.so");
writeControl writeTime;
log true;
writeFields true;
}

Lambda21
{
type Lambda2;
libs ("libfieldFunctionObjects.so");
writeControl writeTime;
log true;
writeFields true;
}
}
  • 基于OF4重新运行求解器(如 interFoam, pisoFoam 等),或运行后处理命令提取已有的结果:
1
2
postProcess -func Q
postProcess -func Lambda2

或者 基于 https://www.cfd-online.com/Forums/openfoam-solving/190302-how-use-postprocess-decomposed-parallel-case.html

1
mpirun -np 2 pimpleFoam -parallel -postProcess .....

2. paraview 后处理

Pasted image 20260715222245.png

第二步:在 ParaView 中显示涡等值面

  • 在终端中进入您的算例目录,运行 paraFoam 命令打开 ParaView。
  • 在左侧面板的 Volume Fields 中勾选 QLambda2(确保读取了这些字段的数据),点击 Apply
  • 在顶部工具栏的过滤菜单中选择 Contour(等值面) 滤镜。
  • 在左侧属性面板的 Contour by 下拉菜单中,选择 QLambda2
  • 设置等值面的数值(Value):
    • 对于 Q 准则:输入一个正数(例如 10 或 100),具体数值取决于您的流动特征,它代表旋转大于剪切应变的区域。
    • 对于 λ₂ 准则:输入一个负数(例如 -0.01 或 -10),负值区域代表压力极小值点,能很好地勾勒出涡核。
  • 点击 Apply 即可生成三维涡结构等值面。可以配合染色(Color by UQ)以增强显示效果。

涡结构提取

转自https://xiaopingqiu.github.io/2016/05/22/QAndLambda/

为了研究湍流的涡结构,需要有一些方法来将涡结构提取出来,比图在文章中常见类似这种图:
涡结构涡结构

本篇介绍怎么在 OpenFOAM 中提取涡结构。

历史上曾用过的涡结构提取有以下几种:

  1. 压强的局部极小值
    在形成涡的地方,通常伴随着压强的极小值。比如:
    img
    这种方法的缺点在于,缺乏客观的压力阈值来捕捉所有的涡结构,而且,压力出现极值的地方不见得就真的有涡。

  2. 流线
    通过流线的封闭来显示涡的结构也是一种常见方法,比如
    img
    这种方法有一个最明显的缺点是,流线不满足伽利略不变性,即,如果换一个参考系,则可能显示出来的“涡结构”就完全不一样了。另外,这种方法也难以分辨两个很靠近的涡。

  3. 涡量的模
    用涡量的模来显示涡结构是一种很常用的方法,类似这样
    img
    这种方法在自由剪切流中很有效,不过,对于壁面束缚流动则不太适用,原因是背景流动的剪切性导致的涡量模可以达到跟涡结构处的涡量的模差不多大小,这就使得涡结构难以从背景流动中分离出来了。并且,涡量的模的最大值通常发生在壁面上,而涡的核心显然不可能出现在壁面上。所以这种方法不适合用于提取边界层附近的涡结构。

    OpenFOAM 中提供了两种方法来提取涡结构:Q 和 Lambda2。

  • 速度梯度张量的二阶不变量
    速度梯度 $\nabla \mathbf{U} $的二阶不变量 Q 的定义为
    $$
    Q = \frac{1}{2}\Big ( ||\mathbf{W}||^2 - ||\mathbf{S}||^2 \Big )
    $$

    其中

    $\mathbf{W} = \frac{1}{2} \Big ( \nabla \mathbf{U} - (\nabla \mathbf{U}) ^{\mathrm{T}} \Big ) \ ||\mathbf{W}|| = (\mathbf{W}:\mathbf{W})^{1/2} \ \mathbf{S} = \frac{1}{2} \Big ( \nabla \mathbf{U} + (\nabla \mathbf{U}) ^{\mathrm{T}} \Big ) \ ||\mathbf{S}|| = (\mathbf{S}:\mathbf{S})^{1/2}$

可以用 Q > 0 来作为涡结构存在的盘踞。
在 OpenFOAM 中,有一个程序用来计算 Q,名字就叫 Q。在流场计算完毕以后,可以运行 Q,然后在 paraview 中显示 Q 值大于 0 的等值面来显示涡的结构。只是,OpenFOAM 中 Q 的计算用的是另一种方法:

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
//Q.C 
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))))
);

代码里注释说这是另一种计算 Q 的方法,与上面公式的计算方法差别很小。

  • 张量 $\mathbf{W} \cdot \mathbf{W} + \mathbf{S} \cdot \mathbf{S}$ 的第二大特征值

另一种判据是 $\mathbf{W} \cdot \mathbf{W} + \mathbf{S} \cdot \mathbf{S}$ 的第二大特征值$ \lambda _ 2 < 0$。
在 OpenFOAM 中有一个程序用来计算 $\lambda _ 2$ :Lambda2

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
//Lambda2.C
const volTensorField gradU(fvc::grad(U));

volTensorField SSplusWW
(
(symm(gradU) & symm(gradU)) + (skew(gradU) & skew(gradU))
);

volScalarField Lambda2
(
IOobject
(
"Lambda2",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::NO_WRITE
),
-eigenValues(SSplusWW)().component(vector::Y)
);

Info<< " Writing -Lambda2" << endl;
Lambda2.write();

注意,OpenFOAM 返回的是 $- \lambda _ 2$,所以,在计算了 Lambda2 后,需要通过 Lambda2 大于 0 的等值面来显示涡结构。本篇开头第一张图片,显示的是圆柱绕流的 Lambda2 = 500 等值面。

参考
Eugene de Villiers, The Potential of Large Eddy Simulation for the Modeling of Wall Bounded Flows, Ph.D Thesis, Imperial College of Science, 2005.