很久之前就把网站搭建好了，但一直没往里写东西。最近被一位神秘的大佬催更，写篇文章吧。
无意中看到了一个之前写的后处理工具，决定就拿这个作第一篇博客。

需求

OpenFOAM不能对计算结果进行周向平均处理。所以就写了个后处理程序，用于对轴向对称流场的三维结果进行周向平均。

解决思路

确定一个周向平均操作需要给定两个量，基点basePoint和旋转向量spanwiseVector。接下来可以通过两种方法实现周向平均过程。

  • 对流场的所有网格遍历，将所有位于同一圆周线上的点加和求平均。可以这么实现：对原始坐标系进行坐标变换，将原来的三维坐标(x,y,z)映射到原点为basePoint高度方向为spanwiseVector的圆柱坐标系（r,z,theta）中；然后在theta方向上求平均，即将所有具有相同（r,z）的节点上的物理量取平均。

  • 取有限个经过基点basePoint和旋转向量spanwiseVector的截面，将各截面上的物理量值求平均。这么实现：过基点和旋转向量选取任意一个截面，提取截面上各点的物理量值；将该初始截面旋转一个角度dTheta，对旋转后的新的截面也提取一次。如此，旋转一周均匀截取n个截面。最终对这n个截面上的物理量值取平均，得到周向平均的结果。

显然，当截面数n比较小的情况下，方案二的计算量显著小于方案一；和理论值之间的偏差取决于截面数量，精度可控。所以接下来就按照方案二的思路来介绍，代码是基础OpenFOAM-4.1写的。

具体步骤

  1. 首先需要确定一个初始截面，提取截面上的坐标，在圆柱坐标系中记为（r,z,0）；
  2. 将初始截面上的每个点绕基点和旋转向量旋转角度dTheta，得到一个旋转后的采样点的坐标（r,z,theta）；
  3. 通过插值的方法得到旋转后的采样点上的物理量值。
  4. 重复步骤2-3，依次获取n个采样点的物理量值，取平均得到初始截面上任意一点（r,z,0）的周向平均值。

实现方法上参考了这篇cfd-online帖子。

1. 确定初始截面

值得一提的是，周向平均并不是一个普适性的后处理过程，它对计算域的几何形状有要求。有些几何形状并不是轴向对称，或者不是关于用户提供的基点basePoint和旋转向量normVector轴向对称的，从原理上讲本就无法进行周向平均。那么这样的几何形状怎么处理，这里涉及到一个容错的问题，这在后面会提到，这应该也是为什么OpenFOAM官方没有提供周向平均后处理程序的一个原因。

本着尽可能用OpenFOAM现有类和函数的原则，我们找到了cuttingPlane类和sampledSurface类，前者可以根据基点basePoint和法向量normVector截取一个初始截面并获得在该截面上的点的坐标，后者则是用来获取这些点的物理量值的。显然，我们直接继承这两个类就可以完成步骤1的主要内容了，而OpenFOAM恰好有一个很好的例子，就是sampledPlan类，同时继承了上面提到的两个类，它在这里

1
$FOAM_SRC/sampling/sampledSurface/sampledPlan

所以我们就按照sampledPlan类改写，新建一个新的类sampledPlaneSpanwise，这里截取一段类的声明，新增了采样点个数nPoints_和旋转向量spanwiseVector_. 基点basePoint和初始截面的法向量normVector就不需要了，因为他们可以直接从cuttingPlane类继承来。

sampledPlaneSpanwise/sampledPlaneSpanwise.H

1
2
3
4
5
6
7
8
9
10
11
12
13
14
class sampledPlaneSpanwise
:
    public sampledSurface,
    public cuttingPlane
{
    // Private data
        //- Direction normal to the plane
            word axis_;

        //- Number of points used to compute the average
            label nPoints_;

        //- Number of points used to compute the average
            vector spanwiseVector_;

在sampledPlaneSpanwise类的构造函数中，有这么一段，用于判断初始截面的法向量normVector和旋转向量spanwiseVector是否垂直，这是容错机制的一部分。

sampledPlaneSpanwise/sampledPlaneSpanwise.C

1
2
3
4
5
6
7
8
9
10
const vector& normal = spanwiseVector_;
const vector& planNorm = this->normal();
if ( mag(normal&planNorm) > SMALL)
{
    FatalErrorIn
        (
             "Foam::sampledAveragePlane::sampledAveragePlane"
        )   << "The spanwiseVector normal is not perpendicular to the plane normal"
            << exit(FatalError);
}

2. 旋转采样

圆周采样类circleSet，可以给定圆心，半径和旋转轴在圆周方向上均匀采样，得到一组采样点的坐标。这么一来，圆柱坐标变换和旋转坐标变换的步骤也省了，所以说还是用OpenFOAM现有类和函数好啊。

circleSet类的声明和定义在这里：

1
$FOAM_SRC/sampling/sampledSet/circle

所以，旋转采样就可以这么实现：为每个初始截面上的点currentPoint，创建一个circleSet对象，其中包含nPoints_个采样点，相邻两个采样点之间的夹角是dTheta。

sampledPlaneSpanwise/sampledPlaneSpanwiseTemplates.C

1
2
// Extract the circle in the homogeneous direction
    circleSet   line("circleLine", mesh(), searchEngine, axis_, circleOrigin, normalVector, currentPoint, dTheta);

circleSet在创建新对象时会调用findCell函数查找改采样点是否在计算域内，不在的话会有输出警告。这也是前面提到的容错机制的一种，只有选取了合适的基点basePoint法向量normVector和旋转向量spanwiseVector才能保证所有的采样点都在计算域内。

3. 插值获取采样点上的物理量

参考 $FOAM_SRC/sampling/sampledSurface/sampledPlan/sampledPlanTemplates.C 当中的实现方法，调用interpolate函数，就可以得到圆周采样点上的物理量值。实现如下

sampledPlaneSpanwise/sampledPlaneSpanwiseTemplates.C

1
2
3
4
5
6
7
8
9
forAll(line, lineI)
{
        linevalues[lineI] = interpolator.interpolate
        (
                line[lineI],
                line.cells()[lineI],
                line.faces()[lineI]
        );
}

2. 取周向平均

得到圆周采样点上的物理量值后，对nPoints_个点求平均即可得到改点处的周向平均值。对初始截面上每个点在每个物理量场(p, U, T …)上作相同地采样和平均处理，即可得到初始截面上的所有点的周向平均值，输出即可。而输出的格式因为继承了sampling类，既可以导出vtk格式又可以得到行列式的raw格式。具体参考sampledPlan的使用方法一样。文末提供的下载链接中也提供了一个算例作参考。

下载链接

我放在github上了https://github.com/ZmengXu/sampledAveragePlane。如果访问不了，可以发邮件找我要。