dsmcFoam solver分析
分析这个solver的主要目的是弄清楚我用这个solver算热驱方腔流动问题为何速度场老是计算得有很大的噪音。即使统计了很长时间。我需要弄清楚这里面的场是如何平均的。
dsmcFoam 的架构：
主要的源文件：
dsmcFoma 源码：
int main(int argc, char *argv[])
{
#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //
Info<< nl << "Constructing dsmcCloud " << endl;
<span class="n">dsmcCloud</span> <span class="n">dsmc</span><span class="p">(</span><span class="s">"dsmc"</span><span class="p">,</span> <span class="n">mesh</span><span class="p">);</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">"</span><span class="se">\n</span><span class="s">Starting time loop</span><span class="se">\n</span><span class="s">"</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="k">while</span> <span class="p">(</span><span class="n">runTime</span><span class="p">.</span><span class="n">loop</span><span class="p">())</span>
<span class="p">{</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">"Time = "</span> <span class="o"><<</span> <span class="n">runTime</span><span class="p">.</span><span class="n">timeName</span><span class="p">()</span> <span class="o"><<</span> <span class="n">nl</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">dsmc</span><span class="p">.</span><span class="n">evolve</span><span class="p">();</span> <span class="c1">//<<------ 这里是演化步
dsmc.info(); //<<------ dump到屏幕一些基本信息 runTime.write(); //<<------- 这里很关键，到底写的是什么 Info<< nl << "ExecutionTime = " << runTime.elapsedCpuTime() << " s" << " ClockTime = " << runTime.elapsedClockTime() << " s" << nl << endl; }
<span class="n">Info</span><span class="o"><<</span> <span class="s">"End</span><span class="se">\n</span><span class="s">"</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="k">return</span><span class="p">(</span><span class="mi">0</span><span class="p">);</span>
}
DsmcCloud 关键函数分析
DsmcCloud::evolve() 主要演化过程在这里
template<class ParcelType>
void Foam::DsmcCloud<ParcelType>::evolve()
{
typename ParcelType::trackingData td(*this);
<span class="c1">// Reset the data collection fields
resetFields(); //<<------ 每个演化步都需要重置 rhoN rhoM linearKE_ ... 等等这些场为0； if (debug) { this->dumpParticlePositions(); }
<span class="c1">// Insert new particles from the inflow boundary
this->inflowBoundary().inflow(); //<<------ 处理边界 // Move the particles ballistically with their current velocities Cloud<ParcelType>::move(td, mesh.time().deltaTValue()); //<<------ 自由迁移，得到particle新的位置 // Update cell occupancy buildCellOccupancy(); //<<------- // Calculate new velocities via stochastic collisions collisions(); //<<------ 碰撞处理，得到particle新的速度 // Calculate the volume field data calculateFields(); //<<------每个演化步都需要统计 rhoN rhoM linearKE ... 等等这些场为； }
DsmcCloud::resetFields() 每个演化步都需要重置这些场为0，为开始统计做准备
template<class ParcelType>
void Foam::DsmcCloud<ParcelType>::resetFields()
{
q_ = dimensionedScalar("zero", dimensionSet(1, 0, -3, 0, 0), 0.0);
<span class="n">fD_</span> <span class="o">=</span> <span class="n">dimensionedVector</span>
<span class="p">(</span>
<span class="s">"zero"</span><span class="p">,</span>
<span class="n">dimensionSet</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">),</span>
<span class="n">vector</span><span class="o">::</span><span class="n">zero</span>
<span class="p">);</span>
<span class="n">rhoN_</span> <span class="o">=</span> <span class="n">dimensionedScalar</span><span class="p">(</span><span class="s">"zero"</span><span class="p">,</span> <span class="n">dimensionSet</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">),</span> <span class="n">VSMALL</span><span class="p">);</span>
<span class="n">rhoM_</span> <span class="o">=</span> <span class="n">dimensionedScalar</span><span class="p">(</span><span class="s">"zero"</span><span class="p">,</span> <span class="n">dimensionSet</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">),</span> <span class="n">VSMALL</span><span class="p">);</span>
<span class="n">dsmcRhoN_</span> <span class="o">=</span> <span class="n">dimensionedScalar</span><span class="p">(</span><span class="s">"zero"</span><span class="p">,</span> <span class="n">dimensionSet</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">),</span> <span class="mf">0.0</span><span class="p">);</span>
<span class="n">linearKE_</span> <span class="o">=</span> <span class="n">dimensionedScalar</span><span class="p">(</span><span class="s">"zero"</span><span class="p">,</span> <span class="n">dimensionSet</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">),</span> <span class="mf">0.0</span><span class="p">);</span>
<span class="n">internalE_</span> <span class="o">=</span> <span class="n">dimensionedScalar</span><span class="p">(</span><span class="s">"zero"</span><span class="p">,</span> <span class="n">dimensionSet</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">),</span> <span class="mf">0.0</span><span class="p">);</span>
<span class="n">iDof_</span> <span class="o">=</span> <span class="n">dimensionedScalar</span><span class="p">(</span><span class="s">"zero"</span><span class="p">,</span> <span class="n">dimensionSet</span><span class="p">(</span><span class="mi">0</span><span class="p">,</span> <span class="o">-</span><span class="mi">3</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">),</span> <span class="n">VSMALL</span><span class="p">);</span>
<span class="n">momentum_</span> <span class="o">=</span> <span class="n">dimensionedVector</span>
<span class="p">(</span>
<span class="s">"zero"</span><span class="p">,</span>
<span class="n">dimensionSet</span><span class="p">(</span><span class="mi">1</span><span class="p">,</span> <span class="o">-</span><span class="mi">2</span><span class="p">,</span> <span class="o">-</span><span class="mi">1</span><span class="p">,</span> <span class="mi">0</span><span class="p">,</span> <span class="mi">0</span><span class="p">),</span>
<span class="n">vector</span><span class="o">::</span><span class="n">zero</span>
<span class="p">);</span>
}
DsmcCloud::calculateFields(): 统计场
template<class ParcelType>
void Foam::DsmcCloud<ParcelType>::calculateFields()
{
scalarField& rhoN = rhoN.internalField();
scalarField& rhoM = rhoM.internalField();
scalarField& dsmcRhoN = dsmcRhoN.internalField();
scalarField& linearKE = linearKE.internalField();
scalarField& internalE = internalE.internalField();
scalarField& iDof = iDof.internalField();
vectorField& momentum = momentum_.internalField();
<span class="n">forAllConstIter</span><span class="p">(</span><span class="k">typename</span> <span class="n">DsmcCloud</span><span class="o"><</span><span class="n">ParcelType</span><span class="o">></span><span class="p">,</span> <span class="o">*</span><span class="k">this</span><span class="p">,</span> <span class="n">iter</span><span class="p">)</span>
<span class="p">{</span>
<span class="k">const</span> <span class="n">ParcelType</span><span class="o">&</span> <span class="n">p</span> <span class="o">=</span> <span class="n">iter</span><span class="p">();</span>
<span class="k">const</span> <span class="n">label</span> <span class="n">cellI</span> <span class="o">=</span> <span class="n">p</span><span class="p">.</span><span class="n">cell</span><span class="p">();</span>
<span class="n">rhoN</span><span class="p">[</span><span class="n">cellI</span><span class="p">]</span><span class="o">++</span><span class="p">;</span>
<span class="n">rhoM</span><span class="p">[</span><span class="n">cellI</span><span class="p">]</span> <span class="o">+=</span> <span class="n">constProps</span><span class="p">(</span><span class="n">p</span><span class="p">.</span><span class="n">typeId</span><span class="p">()).</span><span class="n">mass</span><span class="p">();</span>
<span class="n">dsmcRhoN</span><span class="p">[</span><span class="n">cellI</span><span class="p">]</span><span class="o">++</span><span class="p">;</span>
<span class="n">linearKE</span><span class="p">[</span><span class="n">cellI</span><span class="p">]</span> <span class="o">+=</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">constProps</span><span class="p">(</span><span class="n">p</span><span class="p">.</span><span class="n">typeId</span><span class="p">()).</span><span class="n">mass</span><span class="p">()</span><span class="o">*</span><span class="p">(</span><span class="n">p</span><span class="p">.</span><span class="n">U</span><span class="p">()</span> <span class="o">&</span> <span class="n">p</span><span class="p">.</span><span class="n">U</span><span class="p">());</span>
<span class="n">internalE</span><span class="p">[</span><span class="n">cellI</span><span class="p">]</span> <span class="o">+=</span> <span class="n">p</span><span class="p">.</span><span class="n">Ei</span><span class="p">();</span>
<span class="n">iDof</span><span class="p">[</span><span class="n">cellI</span><span class="p">]</span> <span class="o">+=</span> <span class="n">constProps</span><span class="p">(</span><span class="n">p</span><span class="p">.</span><span class="n">typeId</span><span class="p">()).</span><span class="n">internalDegreesOfFreedom</span><span class="p">();</span>
<span class="n">momentum</span><span class="p">[</span><span class="n">cellI</span><span class="p">]</span> <span class="o">+=</span> <span class="n">constProps</span><span class="p">(</span><span class="n">p</span><span class="p">.</span><span class="n">typeId</span><span class="p">()).</span><span class="n">mass</span><span class="p">()</span><span class="o">*</span><span class="n">p</span><span class="p">.</span><span class="n">U</span><span class="p">();</span>
<span class="p">}</span>
<span class="n">rhoN</span> <span class="o">*=</span> <span class="n">nParticle_</span><span class="o">/</span><span class="n">mesh</span><span class="p">().</span><span class="n">cellVolumes</span><span class="p">();</span>
<span class="n">rhoN_</span><span class="p">.</span><span class="n">correctBoundaryConditions</span><span class="p">();</span>
<span class="n">rhoM</span> <span class="o">*=</span> <span class="n">nParticle_</span><span class="o">/</span><span class="n">mesh</span><span class="p">().</span><span class="n">cellVolumes</span><span class="p">();</span>
<span class="n">rhoM_</span><span class="p">.</span><span class="n">correctBoundaryConditions</span><span class="p">();</span>
<span class="n">dsmcRhoN_</span><span class="p">.</span><span class="n">correctBoundaryConditions</span><span class="p">();</span>
<span class="n">linearKE</span> <span class="o">*=</span> <span class="n">nParticle_</span><span class="o">/</span><span class="n">mesh</span><span class="p">().</span><span class="n">cellVolumes</span><span class="p">();</span>
<span class="n">linearKE_</span><span class="p">.</span><span class="n">correctBoundaryConditions</span><span class="p">();</span>
<span class="n">internalE</span> <span class="o">*=</span> <span class="n">nParticle_</span><span class="o">/</span><span class="n">mesh</span><span class="p">().</span><span class="n">cellVolumes</span><span class="p">();</span>
<span class="n">internalE_</span><span class="p">.</span><span class="n">correctBoundaryConditions</span><span class="p">();</span>
<span class="n">iDof</span> <span class="o">*=</span> <span class="n">nParticle_</span><span class="o">/</span><span class="n">mesh</span><span class="p">().</span><span class="n">cellVolumes</span><span class="p">();</span>
<span class="n">iDof_</span><span class="p">.</span><span class="n">correctBoundaryConditions</span><span class="p">();</span>
<span class="n">momentum</span> <span class="o">*=</span> <span class="n">nParticle_</span><span class="o">/</span><span class="n">mesh</span><span class="p">().</span><span class="n">cellVolumes</span><span class="p">();</span>
<span class="n">momentum_</span><span class="p">.</span><span class="n">correctBoundaryConditions</span><span class="p">();</span>
}
看看DsmcCloud的构造函数，可以发现它主要存储哪些东西：
// * * * * * * * * * * * * * * * * Constructors * * * * * * * * * * * * * * //
template<class ParcelType>
Foam::DsmcCloud<ParcelType>::DsmcCloud
(
const word& cloudName,
const fvMesh& mesh,
bool readFields
)
:
Cloud<ParcelType>(mesh, cloudName, false),
DsmcBaseCloud(),
cloudName(cloudName),
mesh(mesh),
particleProperties
(
IOobject
(
cloudName + "Properties",
mesh.time().constant(),
mesh,
IOobject::MUST_READ_IF_MODIFIED,
IOobject::NO_WRITE
)
),
typeIdList(particleProperties.lookup("typeIdList")),
nParticle(readScalar(particleProperties.lookup("nEquivalentParticles"))),
cellOccupancy(mesh.nCells()),
sigmaTcRMax //<<------ 这个主要是作为计算的时候用，输出并没有意义。
(
IOobject
(
this->name() + "SigmaTcRMax",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ, //<<------ 这里都是MUST_READ
IOobject::AUTO_WRITE
),
mesh
),
collisionSelectionRemainder(mesh.nCells(), 0),
q
(
IOobject
(
"q",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
fD
(
IOobject
(
"fD",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
rhoN
(
IOobject
(
"rhoN",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
rhoM
(
IOobject
(
"rhoM",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
dsmcRhoN
(
IOobject
(
"dsmcRhoN",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
linearKE
(
IOobject
(
"linearKE",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
internalE
(
IOobject
(
"internalE",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
iDof
(
IOobject
(
"iDof",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
momentum
(
IOobject
(
"momentum",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
),
constProps(),
rndGen(label(149382906) + 7183*Pstream::myProcNo()),
boundaryT
(
volScalarField
(
IOobject
(
"boundaryT",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
),
boundaryU
(
volVectorField
(
IOobject
(
"boundaryU",
mesh.time().timeName(),
mesh,
IOobject::MUST_READ,
IOobject::AUTO_WRITE
),
mesh
)
),
binaryCollisionModel
(
BinaryCollisionModel<DsmcCloud<ParcelType> >::New
(
particleProperties,
*this
)
),
wallInteractionModel
(
WallInteractionModel<DsmcCloud<ParcelType> >::New
(
particleProperties,
*this
)
),
inflowBoundaryModel
(
InflowBoundaryModel<DsmcCloud<ParcelType> >::New
(
particleProperties_,
*this
)
)
{
buildConstProps();
<span class="n">buildCellOccupancy</span><span class="p">();</span>
<span class="c1">// Initialise the collision selection remainder to a random value between 0
// and 1. forAll(collisionSelectionRemainder, i) { collisionSelectionRemainder[i] = rndGen_.scalar01(); }
<span class="k">if</span> <span class="p">(</span><span class="n">readFields</span><span class="p">)</span>
<span class="p">{</span>
<span class="n">ParcelType</span><span class="o">::</span><span class="n">readFields</span><span class="p">(</span><span class="o">*</span><span class="k">this</span><span class="p">);</span>
<span class="p">}</span>
}
可以发现DsmcCloud的演化过程十分清晰明了。就是内部记录了一些extensive fields，和 particle的位置、速度。然后每个演化步统计这些 extensive fields，然后执行迁移和碰撞，边界处理。这些extensive fields完全由新时刻的particle速度计算出来。
fieldAverage 分析
在一个dsmcFoam的case里面，controlDict 文件里面要设置在求解过程中输出一些时间平均的场, 这由OpenFOAM本身提供的fieldAverage 和dsmcFOAM提供的dsmcField 这两个Function object来完成。
NOTE: 关于fieldAverage Function object到底是怎样进行时间平均的， 这篇博文有很好的分析。
下面这个例子是一个dsmcFoam case里controlDict文件的functions部分:
functions // 这里面是设置了两个时间平均计算，一个是fieldAverage1，另一个是dsmcFields1
{
fieldAverage1
{
type fieldAverage; //<<------ 这个是OpenFOAM本身提供的。
functionObjectLibs ( "libfieldFunctionObjects.so" );
outputControl outputTime; //<<------ 随着每次dsmcFoam输出AUTO_WRITE属性的
//场（由controlDict里面的writeControl控制）
//的时候都会输出这些时均场。
resetOnOutput false;
<span class="n">fields</span>
<span class="p">(</span>
<span class="n">rhoN</span>
<span class="p">{</span>
<span class="n">mean</span> <span class="n">on</span><span class="p">;</span>
<span class="n">prime2Mean</span> <span class="n">off</span><span class="p">;</span>
<span class="n">base</span> <span class="n">time</span><span class="p">;</span> <span class="c1">// 如果base设为time的话，每个演化步都会计算一次时均
//windows 2; //窗口长度为2，就是相对于没有windows,统计频率会降低一倍。 } rhoM { mean on; prime2Mean off; base time; } ... //省略 fD { mean on; prime2Mean off; base time; } ); } dsmcFields1 { type dsmcFields; functionObjectLibs ( "libutilityFunctionObjects.so" ); enabled true; outputControl outputTime; } }
这里面设置每个需要计算时均值的场的设置。包括mean,prome2Mean,base,window。时均场输出的文件名后面会多一个Mean。 如momentum的时均场为momentumMean。
上面的设置中base为time，发现计算过程中每个演化步(DsmcCloud::evolve())都计算了时均场，log里面每个演化都会出现：
fieldAverage fieldAverage1 output:
Calculating averages
- 如果不设置
window: 每次计算的时均场都是从求解开始到目前所有步的时均场。 - 如果设置了
window:base用来指定作时间平均的基础，是基于时间步数(ITER)还是物理时间(time); window用来作平均的时间段的长度，如果不设定，则求的是从开始到当前时间这个时间段的平均值。window的数值的实际含义依base而定，如果base是ITER，则window=20表示当前步及其前 19 个时间步从 20 个时间步内的平均，而如果base是time,则表示的是 20s 内的平均。
在输出数据文件时（每evolve 100次）还会输出fieldAverage1设置的时均场(writing average fields)：
fieldAverage fieldAverage1 output:
Calculating averages
writing average fields //<------ 表示输出时均场
Calculating dsmcFields.
Calculating UMean field.
Calculating translationalT field.
Calculating internalT field.
Calculating overallT field.
Calculating pressure field.
mag(UMean) max/min : 76.50092418551631 1.282028354585403
translationalT max/min : 436.5786449656766 219.0809500397843
internalT max/min : 0 0
overallT max/min : 436.5786449656766 219.0809500397843
p max/min : 67296.44015303567 42358.60866053814
dsmcFields written.
NOTE: 可以发现在writing average fields的时候dsmcFields function object被触发了，它计算UMean,translationalT,internalT,overallT,pressure等场。
可以定义两个开关分别是resetOnOutput和resetOnRestart, 其意义是：
resetOnRestart的值决定当solver继续运行时，是否要读取最近一个时间步的meanField的值来计算接下来时刻的时均值；>resetOnOutput，顾名思义，是否要在每一次输出到文件以后重置meanField的值。这两个开关的默认值都是false。"
这里可能有个疑问，如果不做这两个时均计算会怎样？会不会影响计算过程？毕竟dsmc的计算过程不依赖于这些宏观量场以及它们的时均场。这个疑问是很自然的。可以通过一个例子分析一下：
如果在controlDict文件里面如下设置
startTime 0;
stopAt endTime;
endTime 1e-8;
deltaT 1e-11;
writeControl timeStep;
writeInterval 100;
表示计算1000步，每100步输出。
1. 如果我们再注释掉functions部分, 也就是不求平均。运行dsmcFoam，也没有问题，不会报错。运行完会产生10个数据文件目录（1000/100=10）。
[lhzhu@ws3 dsmc1]$ ls
0 1e-08 1e-09 2e-09 3e-09 4e-09 5e-09 6e-09 7e-09 8e-09 9e-09 constant log PyFoamHistory system
每个目录里面只有当前步的统计出来的场(非时均)：
[lhzhu@ws3 dsmc1]$ ls 1e-09
boundaryT boundaryU dsmcRhoN dsmcSigmaTcRMax fD iDof internalE lagrangian linearKE momentum q rhoM rhoN uniform
- 但是如果我们如果只注释掉
functions里面的fieldAverage1部分，或者只去掉fieldAverage1里面的任何一个时均场，则会报错。比如去掉fD的时均计算，会有Fatal Error如下：
--> FOAM FATAL ERROR:
request for volVectorField fDMean from objectRegistry region0 failed
这个错误应该是在执行dsmcFields function object 时产生的。下面我们分析dsmcFields。
dsmcFields分析
源文件位置： https://github.com/OpenFOAM/OpenFOAM-2.3.x/tree/master/src/postProcessing/functionObjects/utilities/dsmcFields
dsmcFields.H
Description
Calculate intensive fields:
- UMean
- translationalT
- internalT
- overallT
from averaged extensive fields from a DSMC calculation.
SourceFiles
dsmcFields.C
IOdsmcFields.H
class dsmcFields
{
// Private data
//- Name of this set of dsmcFields objects
word name_;
<span class="k">const</span> <span class="n">objectRegistry</span><span class="o">&</span> <span class="n">obr_</span><span class="p">;</span>
<span class="c1">//- on/off switch
bool active_;
<span class="c1">// Private Member Functions
//- Disallow default bitwise copy construct dsmcFields(const dsmcFields&);
<span class="c1">//- Disallow default bitwise assignment
void operator=(const dsmcFields&);
public:
<span class="c1">//- Runtime type information
TypeName("dsmcFields");
<span class="c1">// Constructors
//- Construct for given objectRegistry and dictionary. // Allow the possibility to load fields from files dsmcFields ( const word& name, const objectRegistry&, const dictionary&, const bool loadFromFiles = false );
<span class="c1">//- Destructor
virtual ~dsmcFields();
<span class="c1">// Member Functions
//- Return name of the set of dsmcFields virtual const word& name() const { return name_; }
<span class="c1">//- Read the dsmcFields data
virtual void read(const dictionary&);
<span class="c1">//- Execute, currently does nothing
virtual void execute();
<span class="c1">//- Execute at the final time-loop, currently does nothing
virtual void end();
<span class="c1">//- Called when time was set at the end of the Time::operator++
virtual void timeSet();
<span class="c1">//- Calculate the dsmcFields and write
virtual void write(); //<------ 最重要的是这个函数 //- Update for changes of mesh virtual void updateMesh(const mapPolyMesh&) {}
<span class="c1">//- Update for changes of mesh
virtual void movePoints(const polyMesh&) {} };
dsmcFields::write()
void Foam::dsmcFields::write()
{
if (active_)
{
word rhoNMeanName = "rhoNMean";
word rhoMMeanName = "rhoMMean";
word momentumMeanName = "momentumMean";
word linearKEMeanName = "linearKEMean";
word internalEMeanName = "internalEMean";
word iDofMeanName = "iDofMean";
word fDMeanName = "fDMean";
<span class="k">const</span> <span class="n">volScalarField</span><span class="o">&</span> <span class="n">rhoNMean</span> <span class="o">=</span> <span class="n">obr_</span><span class="p">.</span><span class="n">lookupObject</span><span class="o"><</span><span class="n">volScalarField</span><span class="o">></span>
<span class="p">(</span>
<span class="n">rhoNMeanName</span>
<span class="p">);</span>
<span class="k">const</span> <span class="n">volScalarField</span><span class="o">&</span> <span class="n">rhoMMean</span> <span class="o">=</span> <span class="n">obr_</span><span class="p">.</span><span class="n">lookupObject</span><span class="o"><</span><span class="n">volScalarField</span><span class="o">></span>
<span class="p">(</span>
<span class="n">rhoMMeanName</span>
<span class="p">);</span>
<span class="k">const</span> <span class="n">volVectorField</span><span class="o">&</span> <span class="n">momentumMean</span> <span class="o">=</span> <span class="n">obr_</span><span class="p">.</span><span class="n">lookupObject</span><span class="o"><</span><span class="n">volVectorField</span><span class="o">></span>
<span class="p">(</span>
<span class="n">momentumMeanName</span>
<span class="p">);</span>
<span class="k">const</span> <span class="n">volScalarField</span><span class="o">&</span> <span class="n">linearKEMean</span> <span class="o">=</span> <span class="n">obr_</span><span class="p">.</span><span class="n">lookupObject</span><span class="o"><</span><span class="n">volScalarField</span><span class="o">></span>
<span class="p">(</span>
<span class="n">linearKEMeanName</span>
<span class="p">);</span>
<span class="k">const</span> <span class="n">volScalarField</span><span class="o">&</span> <span class="n">internalEMean</span> <span class="o">=</span> <span class="n">obr_</span><span class="p">.</span><span class="n">lookupObject</span><span class="o"><</span><span class="n">volScalarField</span><span class="o">></span>
<span class="p">(</span>
<span class="n">internalEMeanName</span>
<span class="p">);</span>
<span class="k">const</span> <span class="n">volScalarField</span><span class="o">&</span> <span class="n">iDofMean</span> <span class="o">=</span> <span class="n">obr_</span><span class="p">.</span><span class="n">lookupObject</span><span class="o"><</span><span class="n">volScalarField</span><span class="o">></span>
<span class="p">(</span>
<span class="n">iDofMeanName</span>
<span class="p">);</span>
<span class="k">const</span> <span class="n">volVectorField</span><span class="o">&</span> <span class="n">fDMean</span> <span class="o">=</span> <span class="n">obr_</span><span class="p">.</span><span class="n">lookupObject</span><span class="o"><</span><span class="n">volVectorField</span><span class="o">></span>
<span class="p">(</span>
<span class="n">fDMeanName</span>
<span class="p">);</span>
<span class="k">if</span> <span class="p">(</span><span class="n">min</span><span class="p">(</span><span class="n">mag</span><span class="p">(</span><span class="n">rhoNMean</span><span class="p">)).</span><span class="n">value</span><span class="p">()</span> <span class="o">></span> <span class="n">VSMALL</span><span class="p">)</span>
<span class="p">{</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">"Calculating dsmcFields."</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" Calculating UMean field."</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">volVectorField</span> <span class="n">UMean</span>
<span class="p">(</span>
<span class="n">IOobject</span>
<span class="p">(</span>
<span class="s">"UMean"</span><span class="p">,</span>
<span class="n">obr_</span><span class="p">.</span><span class="n">time</span><span class="p">().</span><span class="n">timeName</span><span class="p">(),</span>
<span class="n">obr_</span><span class="p">,</span>
<span class="n">IOobject</span><span class="o">::</span><span class="n">NO_READ</span>
<span class="p">),</span>
<span class="n">momentumMean</span><span class="o">/</span><span class="n">rhoMMean</span>
<span class="p">);</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" Calculating translationalT field."</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">volScalarField</span> <span class="n">translationalT</span>
<span class="p">(</span>
<span class="n">IOobject</span>
<span class="p">(</span>
<span class="s">"translationalT"</span><span class="p">,</span>
<span class="n">obr_</span><span class="p">.</span><span class="n">time</span><span class="p">().</span><span class="n">timeName</span><span class="p">(),</span>
<span class="n">obr_</span><span class="p">,</span>
<span class="n">IOobject</span><span class="o">::</span><span class="n">NO_READ</span>
<span class="p">),</span>
<span class="mf">2.0</span><span class="o">/</span><span class="p">(</span><span class="mf">3.0</span><span class="o">*</span><span class="n">physicoChemical</span><span class="o">::</span><span class="n">k</span><span class="p">.</span><span class="n">value</span><span class="p">()</span><span class="o">*</span><span class="n">rhoNMean</span><span class="p">)</span>
<span class="o">*</span><span class="p">(</span><span class="n">linearKEMean</span> <span class="o">-</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">rhoMMean</span><span class="o">*</span><span class="p">(</span><span class="n">UMean</span> <span class="o">&</span> <span class="n">UMean</span><span class="p">))</span>
<span class="p">);</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" Calculating internalT field."</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">volScalarField</span> <span class="n">internalT</span>
<span class="p">(</span>
<span class="n">IOobject</span>
<span class="p">(</span>
<span class="s">"internalT"</span><span class="p">,</span>
<span class="n">obr_</span><span class="p">.</span><span class="n">time</span><span class="p">().</span><span class="n">timeName</span><span class="p">(),</span>
<span class="n">obr_</span><span class="p">,</span>
<span class="n">IOobject</span><span class="o">::</span><span class="n">NO_READ</span>
<span class="p">),</span>
<span class="p">(</span><span class="mf">2.0</span><span class="o">/</span><span class="n">physicoChemical</span><span class="o">::</span><span class="n">k</span><span class="p">.</span><span class="n">value</span><span class="p">())</span><span class="o">*</span><span class="p">(</span><span class="n">internalEMean</span><span class="o">/</span><span class="n">iDofMean</span><span class="p">)</span>
<span class="p">);</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" Calculating overallT field."</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">volScalarField</span> <span class="n">overallT</span>
<span class="p">(</span>
<span class="n">IOobject</span>
<span class="p">(</span>
<span class="s">"overallT"</span><span class="p">,</span>
<span class="n">obr_</span><span class="p">.</span><span class="n">time</span><span class="p">().</span><span class="n">timeName</span><span class="p">(),</span>
<span class="n">obr_</span><span class="p">,</span>
<span class="n">IOobject</span><span class="o">::</span><span class="n">NO_READ</span>
<span class="p">),</span>
<span class="mf">2.0</span><span class="o">/</span><span class="p">(</span><span class="n">physicoChemical</span><span class="o">::</span><span class="n">k</span><span class="p">.</span><span class="n">value</span><span class="p">()</span><span class="o">*</span><span class="p">(</span><span class="mf">3.0</span><span class="o">*</span><span class="n">rhoNMean</span> <span class="o">+</span> <span class="n">iDofMean</span><span class="p">))</span>
<span class="o">*</span><span class="p">(</span><span class="n">linearKEMean</span> <span class="o">-</span> <span class="mf">0.5</span><span class="o">*</span><span class="n">rhoMMean</span><span class="o">*</span><span class="p">(</span><span class="n">UMean</span> <span class="o">&</span> <span class="n">UMean</span><span class="p">)</span> <span class="o">+</span> <span class="n">internalEMean</span><span class="p">)</span>
<span class="p">);</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" Calculating pressure field."</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">volScalarField</span> <span class="n">p</span>
<span class="p">(</span>
<span class="n">IOobject</span>
<span class="p">(</span>
<span class="s">"p"</span><span class="p">,</span>
<span class="n">obr_</span><span class="p">.</span><span class="n">time</span><span class="p">().</span><span class="n">timeName</span><span class="p">(),</span>
<span class="n">obr_</span><span class="p">,</span>
<span class="n">IOobject</span><span class="o">::</span><span class="n">NO_READ</span>
<span class="p">),</span>
<span class="n">physicoChemical</span><span class="o">::</span><span class="n">k</span><span class="p">.</span><span class="n">value</span><span class="p">()</span><span class="o">*</span><span class="n">rhoNMean</span><span class="o">*</span><span class="n">translationalT</span>
<span class="p">);</span>
<span class="k">const</span> <span class="n">fvMesh</span><span class="o">&</span> <span class="n">mesh</span> <span class="o">=</span> <span class="n">fDMean</span><span class="p">.</span><span class="n">mesh</span><span class="p">();</span>
<span class="n">forAll</span><span class="p">(</span><span class="n">mesh</span><span class="p">.</span><span class="n">boundaryMesh</span><span class="p">(),</span> <span class="n">i</span><span class="p">)</span>
<span class="p">{</span>
<span class="k">const</span> <span class="n">polyPatch</span><span class="o">&</span> <span class="n">patch</span> <span class="o">=</span> <span class="n">mesh</span><span class="p">.</span><span class="n">boundaryMesh</span><span class="p">()[</span><span class="n">i</span><span class="p">];</span>
<span class="k">if</span> <span class="p">(</span><span class="n">isA</span><span class="o"><</span><span class="n">wallPolyPatch</span><span class="o">></span><span class="p">(</span><span class="n">patch</span><span class="p">))</span>
<span class="p">{</span>
<span class="n">p</span><span class="p">.</span><span class="n">boundaryField</span><span class="p">()[</span><span class="n">i</span><span class="p">]</span> <span class="o">=</span>
<span class="n">fDMean</span><span class="p">.</span><span class="n">boundaryField</span><span class="p">()[</span><span class="n">i</span><span class="p">]</span>
<span class="o">&</span> <span class="p">(</span><span class="n">patch</span><span class="p">.</span><span class="n">faceAreas</span><span class="p">()</span><span class="o">/</span><span class="n">mag</span><span class="p">(</span><span class="n">patch</span><span class="p">.</span><span class="n">faceAreas</span><span class="p">()));</span>
<span class="p">}</span>
<span class="p">}</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" mag(UMean) max/min : "</span>
<span class="o"><<</span> <span class="n">max</span><span class="p">(</span><span class="n">mag</span><span class="p">(</span><span class="n">UMean</span><span class="p">)).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="s">" "</span>
<span class="o"><<</span> <span class="n">min</span><span class="p">(</span><span class="n">mag</span><span class="p">(</span><span class="n">UMean</span><span class="p">)).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" translationalT max/min : "</span>
<span class="o"><<</span> <span class="n">max</span><span class="p">(</span><span class="n">translationalT</span><span class="p">).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="s">" "</span>
<span class="o"><<</span> <span class="n">min</span><span class="p">(</span><span class="n">translationalT</span><span class="p">).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" internalT max/min : "</span>
<span class="o"><<</span> <span class="n">max</span><span class="p">(</span><span class="n">internalT</span><span class="p">).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="s">" "</span>
<span class="o"><<</span> <span class="n">min</span><span class="p">(</span><span class="n">internalT</span><span class="p">).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" overallT max/min : "</span>
<span class="o"><<</span> <span class="n">max</span><span class="p">(</span><span class="n">overallT</span><span class="p">).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="s">" "</span>
<span class="o"><<</span> <span class="n">min</span><span class="p">(</span><span class="n">overallT</span><span class="p">).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">" p max/min : "</span>
<span class="o"><<</span> <span class="n">max</span><span class="p">(</span><span class="n">p</span><span class="p">).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="s">" "</span>
<span class="o"><<</span> <span class="n">min</span><span class="p">(</span><span class="n">p</span><span class="p">).</span><span class="n">value</span><span class="p">()</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="n">UMean</span><span class="p">.</span><span class="n">write</span><span class="p">();</span>
<span class="n">translationalT</span><span class="p">.</span><span class="n">write</span><span class="p">();</span>
<span class="n">internalT</span><span class="p">.</span><span class="n">write</span><span class="p">();</span>
<span class="n">overallT</span><span class="p">.</span><span class="n">write</span><span class="p">();</span>
<span class="n">p</span><span class="p">.</span><span class="n">write</span><span class="p">();</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">"dsmcFields written."</span> <span class="o"><<</span> <span class="n">nl</span> <span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="p">}</span>
<span class="k">else</span>
<span class="p">{</span>
<span class="n">Info</span><span class="o"><<</span> <span class="s">"Small value ("</span> <span class="o"><<</span> <span class="n">min</span><span class="p">(</span><span class="n">mag</span><span class="p">(</span><span class="n">rhoNMean</span><span class="p">))</span>
<span class="o"><<</span> <span class="s">") found in rhoNMean field. "</span>
<span class="o"><<</span> <span class="s">"Not calculating dsmcFields to avoid division by zero."</span>
<span class="o"><<</span> <span class="n">endl</span><span class="p">;</span>
<span class="p">}</span>
<span class="p">}</span>
}
明显可以发现dsmcFields所做的只是根据fieldAverage1输出的时均场来求出其他时均场，例如根据时均密度场和时均动量场来求出时均速度场。这个工作也可以由后处理工具dsmcFieldsCalc完成，而不是写成controlDict里面的function object。
总结
综上分析，有以下总结：
- 用
fieldAverage对宏观量场求时均是dsmcFoam里面输出宏观量场的最重要步骤。 -
window控制对那一段时间求平均。如果不设置window就是从开始到当前所有步的平均。 如果设置window, 结合base(ITER/time)，可以设置 从当前时间/步 往前多少s/步 求时均值。 -
dsmcFields辅助计算其他时均场如UMean,overallT等等，它只是对fieldAverage输出的时均场做简单的加减乘除计算。
如何解决我的问题？
要算低速问题(Ma < 0.001)，window必须设置得足够大, 统计得出的速度场噪音才会足够小。
比如我算一个二维问题，计算区域为一个边长为1e-6m的方腔。上壁面温度为200K, 其他三个壁面温度为400K，壁面为漫反射边界。我的设置如下：
等我算对了，再回来更新吧！
