#201923
接下来的一段时间我们将重点讲 div 离散格式。对于离散格式的指定，在 Code of the Week #5: div(phi,U) Gauss linear 中已经进行了说明，此处对相应离散格式进行分析。OpenFOAM 中， \nabla\cdot(UU) 表示为
phi 指的是面上的通量，可以是体积通量，也可以是质量通量，对于体积通量，
对于质量通量，是
div(phi,U) 的向量表达式为 \nabla\cdot(U\otimes U) ，经高斯定理，
因此，div(phi,U) 中的 phi 表示 UdS 网格面的通量，在 OpenFOAM 中是已知量， div(phi,U) 离散格式处理的是 U，从而形成以 U_x (或 U_y, U_z) 为未知量的矩阵。
在 Code of the Week #1: div(phi,U) 中对 linear 离散格式进行了说明。但没有进行详细的解析，现在根据数值计算方面的要求，进行必要的考量和分析。
首先，我们看一下 OpenFOAM 对这一量是如何声明的？
fvVectorMatrix UEqn
(
fvm::ddt(U)
+ fvm::div(phi,U)
- fvm::laplacian(nu,U)
);
每一项都返回 fvVectorMatrix 类型。因此，可以对 fvm::div(phi,U) 进行单独分析。fvm::div() 对应的类定义 (声明/构造) 为 (位于src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.H)
// Constructors
//- Construct given a field to solve for
fvMatrix
(
GeometricField<Type, fvPatchField, volMesh>&,
const dimensionSet&
);
//- Construct as copy
fvMatrix(const fvMatrix<Type>&);
//- Construct as copy of tmp<fvMatrix<Type> > deleting argument
# ifdef ConstructFromTmp
fvMatrix(const tmp<fvMatrix<Type> >&);
# endif
//- Construct from Istream given field to solve for
fvMatrix(GeometricField<Type, fvPatchField, volMesh>&, Istream&);
//- Clone
tmp<fvMatrix<Type> > clone() const
{
return tmp<fvMatrix<Type> >
(
new fvMatrix<Type>(*this)
);
}
可以看出， fvMatrix 可以由4种方式构造：
- GeometryField 和 dimensionedSet
fvMatrix<Type>- 临时的
fmMatrix<Type> - GeometryField 和 Istream 类
fvm::div() 则用第二种和第三种构造函数。fvm::div()函数，fvm 表示的是名称空间，::代表的是作用域，意思是fvm::div中的 div 来自于 fvm 空间。
fvm::div(phi,U) 中有两个参数，其中一个是 phi (surfaceScalarField 类)，另一个是 U (fvVectorField 类)，因此，该函数则调用
div(surfaceScalarField& phi, fvVectorField& vf)
因此，查找 fvm 名称空间，有
template<class Type >
tmp< fvMatrix< Type > > div (const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
template<class Type >
tmp< fvMatrix< Type > > div (const tmp< surfaceScalarField > &tflux, const GeometricField< Type, fvPatchField, volMesh > &vf, const word &name)
template<class Type >
tmp< fvMatrix< Type > > div (const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf)
可以找到三个函数，其中有两个参数的是
template<class Type >
tmp< fvMatrix< Type > > div (const surfaceScalarField &flux, const GeometricField< Type, fvPatchField, volMesh > &vf)
构造函数为 (fvmDiv.C)
template<class Type>
tmp<fvMatrix<Type> >
div
(
const surfaceScalarField& flux,
GeometricField<Type, fvPatchField, volMesh>& vf
)
{
return fvm::div(flux, vf, "div("+flux.name()+','+vf.name()+')');
}
可以看到，这个构造函数也调用了另一个函数：
template<class Type>
tmp<fvMatrix<Type> >
div
(
const surfaceScalarField& flux,
GeometricField<Type, fvPatchField, volMesh>& vf,
const word& name
)
{
return fv::convectionScheme<Type>::New
(
vf.mesh(),
flux,
vf.mesh().schemesDict().divScheme(name) // 第 3 参数，指定格式
)().fvmDiv(flux, vf);
}
这个函数则调用了 New()，而它属于 fv::convectionScheme<Type> 类。而生成的对象调用子函数 fvmDiv()。
首先，看一下 New() 函数 (src/finiteVolume/finiteVolume/convectionSchemes/convectionScheme /convectionScheme.H)
//- Return a pointer to a new convectionScheme created on freestore
static tmp<convectionScheme<Type> > New
(
const fvMesh& mesh,
const surfaceScalarField& faceFlux,
Istream& schemeData
);
类本身还有虚函数
virtual tmp<fvMatrix<Type> > fvmDiv
(
const surfaceScalarField&,
GeometricField<Type, fvPatchField, volMesh>&
) const = 0;
virtual tmp<GeometricField<Type, fvPatchField, volMesh> > fvcDiv
(
const surfaceScalarField&,
const GeometricField<Type, fvPatchField, volMesh>&
) const = 0;
而使用的 New() 函数的代码为
template<class Type>
tmp<convectionScheme<Type> > convectionScheme<Type>::New
(
const fvMesh& mesh,
const surfaceScalarField& faceFlux,
Istream& schemeData
)
{
if (fv::debug)
{
Info<< "convectionScheme<Type>::New"
"(const fvMesh&, const surfaceScalarField&, Istream&) : "
"constructing convectionScheme<Type>"
<< endl;
}
if (schemeData.eof())
{
FatalIOErrorIn
(
"convectionScheme<Type>::New"
"(const fvMesh&, const surfaceScalarField&, Istream&)",
schemeData
) << "Convection scheme not specified" << endl << endl
<< "Valid convection schemes are :" << endl
<< IstreamConstructorTablePtr_->sortedToc()
<< exit(FatalIOError);
}
word schemeName(schemeData);
typename IstreamConstructorTable::iterator cstrIter =
IstreamConstructorTablePtr_->find(schemeName);
if (cstrIter == IstreamConstructorTablePtr_->end())
{
FatalIOErrorIn
(
"convectionScheme<Type>::New"
"(const fvMesh&, const surfaceScalarField&, Istream&)",
schemeData
) << "unknown convection scheme " << schemeName << endl << endl
<< "Valid convection schemes are :" << endl
<< IstreamConstructorTablePtr_->sortedToc()
<< exit(FatalIOError);
}
return cstrIter()(mesh, faceFlux, schemeData);
}
其中的第三个参数 schemeData 来自 vf.mesh().divScheme(name)，实际上， mesh (fvMesh 类) 对象并没有直接的 divScheme 函数，而是继承了 fvSchemes 类，其中包括了 divScheme() (src/finiteVolume/finiteVolume/fvSchemes/fvSchemes.C)：
Foam::ITstream& Foam::fvSchemes::divScheme(const word& name) const
{
if (debug)
{
Info<< "Lookup divScheme for " << name << endl;
}
if (divSchemes_.found(name) || defaultDivScheme_.empty())
{
return divSchemes_.lookup(name);
}
else
{
const_cast<ITstream&>(defaultDivScheme_).rewind();
return const_cast<ITstream&>(defaultDivScheme_);
}
}
其意思是：如果在 system/fvSchemes 文件中找到 div(phi,U) 或 default 项为 none，则返回 div(phi,U) 指定的格式，否则按 default 指定的格式。可以进一步找 lookup() 函数。此处就不展开了。
返回一个 convectionScheme 对象 cstrIter，可能的方式有
- boundedConvectionScheme
- gaussConvectionScheme
- multivariateGaussConvectionScheme
当字典文件写的是 Gauss 时，返回的对象是 gaussConvectionScheme，以 fvmDiv 为例，代码为
tmp<fvMatrix<Type> >
gaussConvectionScheme<Type>::fvmDiv
(
const surfaceScalarField& faceFlux,
GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
tmp<surfaceScalarField> tweights = tinterpScheme_().weights(vf);
const surfaceScalarField& weights = tweights();
tmp<fvMatrix<Type> > tfvm
(
new fvMatrix<Type>
(
vf,
faceFlux.dimensions()*vf.dimensions()
)
);
fvMatrix<Type>& fvm = tfvm();
fvm.lower() = -weights.internalField()*faceFlux.internalField();
fvm.upper() = fvm.lower() + faceFlux.internalField();
fvm.negSumDiag();
forAll(fvm.psi().boundaryField(), patchI)
{
const fvPatchField<Type>& psf = fvm.psi().boundaryField()[patchI];
const fvsPatchScalarField& patchFlux = faceFlux.boundaryField()[patchI];
const fvsPatchScalarField& pw = weights.boundaryField()[patchI];
fvm.internalCoeffs()[patchI] = patchFlux*psf.valueInternalCoeffs(pw);
fvm.boundaryCoeffs()[patchI] = -patchFlux*psf.valueBoundaryCoeffs(pw);
}
if (tinterpScheme_().corrected())
{
fvm += fvc::surfaceIntegrate(faceFlux*tinterpScheme_().correction(vf));
}
return tfvm;
}
这个函数有一个类成员 tinterpScheme_ 指定了面上的量的插值方式 system/fvSchemes。
假设此处使用的是 linear 函数，怎样去打捞这个函数呢？在这里 src/finiteVolume/interpolation/surfaceInterpolation，而 linear 则位于src/finiteVolume/interpolation/surfaceInterpolation/schemes/linear，从上一段可以看出，fvmDiv() 使用了以下几个参数
- tinterpScheme_().weights(vf)
- tinterpScheme_().corrected()
- tinterpScheme_().correction()
weights() 使用一个参数的成员函数，在 linear.H 中其定义为
//- Return the interpolation weighting factors
tmp<surfaceScalarField> weights
(
const GeometricField<Type, fvPatchField, volMesh>&
) const
{
return this->mesh().surfaceInterpolation::weights();
}
weights() 又指向了 makeWeights() (surfaceInterpolation.C)：
const Foam::surfaceScalarField& Foam::surfaceInterpolation::weights() const
{
if (!weightingFactors_)
{
makeWeights();
}
return (*weightingFactors_);
}
makeWeights() 在 surfaceInterpolation.C 定义：
void Foam::surfaceInterpolation::makeWeights() const
{
if (debug)
{
Info<< "surfaceInterpolation::makeWeights() : "
<< "Constructing weighting factors for face interpolation"
<< endl;
}
weightingFactors_ = new surfaceScalarField
(
IOobject
(
"weightingFactors",
mesh_.pointsInstance(),
mesh_
),
mesh_,
dimless
);
surfaceScalarField& weightingFactors = *weightingFactors_;
// Set local references to mesh data
// (note that we should not use fvMesh sliced fields at this point yet
// since this causes a loop when generating weighting factors in
// coupledFvPatchField evaluation phase)
const unallocLabelList& owner = mesh_.owner();
const unallocLabelList& neighbour = mesh_.neighbour();
const vectorField& Cf = mesh_.faceCentres(); // 面的中心
const vectorField& C = mesh_.cellCentres(); // 网格中心
const vectorField& Sf = mesh_.faceAreas(); // 面的面积
// ... and reference to the internal field of the weighting factors
scalarField& w = weightingFactors.internalField();
forAll (owner, facei) // 关键部分
{
// Note: mag in the dot-product.
// For all valid meshes, the non-orthogonality will be less than
// 90 deg and the dot-product will be positive. For invalid
// meshes (d & s <= 0), this will stabilise the calculation
// but the result will be poor.
scalar SfdOwn = mag(Sf[facei] & (Cf[facei] - C[owner[facei]])); // 本网格 (共享) 面到中心的距离 A
scalar SfdNei = mag(Sf[facei] & (C[neighbour[facei]] - Cf[facei])); // 相邻网格中心到共享面的距离 B
w[facei] = SfdNei/(SfdOwn + SfdNei); // B/(A+B)
}
forAll (mesh_.boundary(), patchi) // 边界
{
mesh_.boundary()[patchi].makeWeights
(
weightingFactors.boundaryField()[patchi]
);
}
if (debug)
{
Info<< "surfaceInterpolation::makeWeights() : "
<< "Finished constructing weighting factors for face interpolation"
<< endl;
}
}
从上面的代码可见，内部和边界分别生成 weights 值。内部 weights 值跟网格的面到中心的距离有关，这个应该很好理解吧：
而边界上的值则单独处理，其源程序位于 fvPatch.C：
void Foam::fvPatch::makeWeights(scalarField& w) const
{
w = 1.0;
}
而 linear 插值方式没有 correction 方式。