OpenFOAM

foamTimeAverage

在前面的一篇博文中，我介绍了fieldAverage这个functionObject的用法，其中提到， 可以用window这个参数来控制所计算的时均值的时间范围。如果 base = time且window = 10，那从第10s以后，每个时刻 t 输出的时均值其实相当于从 t-10 到 t 这个时间段内的时均值。但是，根据分析可以发现，这个时均值并不严格等价于从 t-10 时刻到 t 时刻某个场的时均值。有时候，需要从某个时刻才开始计算时均值，而fieldAverage没有参数可以控制从某个时刻才开始计算时间平均。于是我参照OpenFOAM的patchAverage.C的代码写了一个后处理程序，用来计算指定时间段内的某个场的时均值。

foamTimeAverage 简介

foamTimeAverage 是一段简单的后处理程序，其功能是在算例运行结束以后，根据指定的时间段，从数据文件夹里循环读入指定时间段内指定场的数据，并计算该段时间的该指定场的时均值，结果将会输出到该时间段最后一个时刻的数据文件夹内。举例说：

1
foamTimeAverage p -time 0.4:0.5

将会计算0.4 s-0.5 s时间段内p的时均值，结果将输出到0.5内，时均值文件名为p_mean。这个程序只支持计算volField的时均值，不支持surfaceField。目前，这个程序只在OpenFOAM-2.1.1和OpenFOAM-2.3.1下通过可编译测试，运行的结果目前只在OpenFOAM-2.1.1下对volScalarField和volVectorField进行了测试。理论上讲，应该在其他版本的OpenFOAM下也可以正常编译和运行。

以下是foamTimeAverage的代码：

1
2
3
4
5
6
7
8
9
10
11
12
13
14
15
16
17
18
19
20
21
22
23
24
25
26
27
28
29
30
31
32
33
34
35
36
37
38
39
40
41
42
43
44
45
46
47
48
49
50
51
52
53
54
55
56
57
58
59
60
61
62
63
64
65
66
67
68
69
70
71
72
73
74
75
76
77
78
79
80
81
82
83
84
85
86
87
88
89
90
91
92
93
94
95
96
97
98
99
100
101
102
103
104
105
106
107
108
109
110
111
112
113
114
115
116
117
118
119
120
121
122
123
124
125
126
127
128
129
130
131
132
133
134
135
136
137
138
139
140
141
142
143
144
145
146
147
148
149
150
151
152
153
154
155
156
157
158
159
160
161
162
163
164
165
166
167
168
169
170
171
172
173
174
/*---------------------------------------------------------------------------*\
 *=========                 |
 *\\      /  F ield         | OpenFOAM: The Open Source CFD Toolbox
 *\\     /   O peration     |
 *\\    /    A nd           | Copyright (C) 1991-2009 OpenCFD Ltd.
 *\\   /     M anipulation  |
 * -------------------------------------------------------------------------------
 license
 This file is part of OpenFOAM.

 OpenFOAM is free software; you can redistribute it and/or modify it
 under the terms of the GNU General Public License as published by the
 Free Software Foundation; either version 2 of the License, or (at your
 option) any later version.

 OpenFOAM is distributed in the hope that it will be useful, but WITHOUT
 ANY WARRANTY; without even the implied warranty of MERCHANTABILITY or
 FITNESS FOR A PARTICULAR PURPOSE.  See the GNU General Public License
 for more details.

 You should have received a copy of the GNU General Public License
 along with OpenFOAM; if not, write to the Free Software Foundation,
 Inc., 51 Franklin St, Fifth Floor, Boston, MA 02110-1301 USA

Application
   foamTimeAverage

Author
  Xiaoping Qiu 
  q.giskard@gmail.com

Description
Calculates the time average  of the specified volField over the specified time range.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

template<class FieldType>
void calcTimeAverage
(
    fvMesh& mesh,
    const IOobject& fieldHeader,
    const word& fieldName,
    Time& runTime,
    instantList& timeDirs,
    bool& done
)
{
    label nfield = 0;
    const word meanFieldName = fieldName + "_mean"; 

    //Info << "class name = " << fieldHeader.headerClassName() << endl;
    //Info << "typeName = " << FieldType::typeName << endl;
    if(!done && fieldHeader.headerClassName() == FieldType::typeName)
    {
	FieldType dummy
	    (
		IOobject
		(
		    fieldName,
		    runTime.timeName(),
		    mesh,
		    IOobject::MUST_READ
		),
		mesh
	    );

	FieldType meanField
	    (
		IOobject
		(
		    meanFieldName,
		    runTime.timeName(),
		    mesh,
		    IOobject::NO_READ
		),
		dummy
	    );

	meanField *= scalar(0.0);

	forAll(timeDirs, timeI)
	{
	    runTime.setTime(timeDirs[timeI], timeI);
	    Info << "Time = " << runTime.timeName() <<endl;

	    IOobject io
		(
		    fieldName,
		    runTime.timeName(),
		    mesh,
		    IOobject::MUST_READ
		);

	    if (io.headerOk())
	    {
		mesh.readUpdate();

		if(!done && io.headerClassName() == FieldType::typeName)
		{
		    Info << "   Reading " << io.headerClassName() << " " <<io.name() << endl;

		    FieldType field(io, mesh);

		    meanField += field;
		    nfield++;
		}
	    }
	    else
	    {
		Info << "   No Field " << fieldName << endl; 
	    }
	}

	if(nfield > 0)
	{
	    Info << "number of field = " << nfield << endl;
	    meanField /= nfield;
	}

	Info<< "writing to timeDir " << runTime.timeName()  << endl;
	meanField.write();
	done = true;

    }
}
// Main program:

int main(int argc, char *argv[])
{

    Foam::timeSelector::addOptions();
    #include "addRegionOption.H"
    Foam::argList::validArgs.append("fieldName");

#   include "setRootCase.H"
#   include "createTime.H"
    instantList timeDirs = timeSelector::select0(runTime, args);
    runTime.setTime(timeDirs[0], 0);
#   include "createNamedMesh.H"

    // get filename from command line
    const word fieldName = args[1];
    bool done = false;

    IOobject fieldHeader
	(
	    fieldName,
	    runTime.timeName(),
	    mesh,
	    IOobject::MUST_READ
	);
    if(fieldHeader.headerOk())//very important!
    {
	calcTimeAverage<volScalarField>(mesh, fieldHeader, fieldName, runTime, timeDirs, done);
	calcTimeAverage<volVectorField>(mesh, fieldHeader, fieldName, runTime, timeDirs, done);
	calcTimeAverage<volTensorField>(mesh, fieldHeader, fieldName, runTime, timeDirs, done);
	calcTimeAverage<volSymmTensorField>(mesh, fieldHeader, fieldName, runTime, timeDirs, done);
	calcTimeAverage<volSphericalTensorField>(mesh, fieldHeader, fieldName, runTime, timeDirs, done);
    }
    else
    {
	Info<< " Error! No field " << fieldName << endl;
    }

    Info<< "End\n" << endl;

    return 0;
}

// ************************************************************************* //

更详细的参考，包括编译和使用方法，见代码的github页面。

一点说明

foamTimeAverage 这个程序有一个缺点，那就是只能根据输出的文件来计算时均值。举例说，假如算例的时间步长取0.0001 s，每0.01 s 输出一次，那0.4 s-0.5 s这个时间段的时均值就只涉及到0.40s, 0.41s, 0.42s, ... 0.50s这些时刻的值。有时候，如果想计算的是更精确的时均值，即0.4000s, 0.4001s, 0.4002s ... 0.5000s这些时刻的时均值，该怎么办呢？我没有实际测试过，但有一种方法我觉得可行，需要结合fieldAverage和foamTimeAverage来实现，具体操作如下：

  1. 设置fieldAverage，注意将resetOnOutput设置为true，outputControl设置为outputTime，base设为time，window不需要设置，这样，应该就会在每一个输出的时间文件夹里输出一个时均值，这个时均值计算的上一次输出到下一次输出之间的时间段内的时均值，用上面的例子来说，即0.5文件夹内输出的p的时均值pMean是0.4001s, 0.4002s, ... 0.5000s这些时刻的p的时均值。

  2. 用foamTimeAverage计算指定时间段内pMean的时均值。举例说，
    foamTimeAverage pMean -time 0.04:0.05 计算的将是 0.3001s, 0.3002s, ... 0.5000s这些时刻的时均值，这是因为，0.04 文件夹内的pMean是 0.3001s, 0.3002s, ... 0.4000s的时均值，0.05 文件夹内的pMean是0.4001s, 0.4002s, ... 0.5000s时刻内的时均值。

    用以上方法就可以实现计算指定时间段内同时精度更高的时均值了。再次申明一下，上述方法仅仅是我根据原理进行的推演，没有经过检验。