Green-Gauss法梯度计算
计算梯度的核心函数位于/opt/openfoam9/src/finiteVolume/finiteVolume/gradSchemes/gradScheme
目录下的gaussGrad.C文件中的gradf函数
。根据笔者的理解,添加注释后,将代码列于下方:
template<class Type>//Type可以是任何类型,如果是scalar那么就是求标量场的梯度,如果是vector那么就是求矢量场梯度(结果是张量场)
Foam::tmp//函数输出tmp,为了方便C++内存管理,在openFoam中tmp代表任何指定的对象
<
Foam::GeometricField//openFoam中最经常使用的场数据,派生自Field,不仅包含场数值,还包含网格数据、单位和边界条件
<
typename Foam::outerProduct<Foam::vector, Type>::type,//这里输出场的类型不是指定的,而是推导得到的。使用外积的概念可知,如果对一个标量场求梯度,那么会得到矢量场。如果对矢量场求梯度,那么会得到张量场,以此类推...
Foam::fvPatchField,
Foam::volMesh//由fvPatchField和volMesh可知,输出场的数值是存储在网格体心的,而不是面心或者顶点
>
>
Foam::fv::gaussGrad<Type>::gradf
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf,//gradf函数的第一个参数:是一个待求梯度的场,注意:由fvsPatchField和surfaceMesh可知,该场的数值是存储在网格面心的,而不是体心或者顶点。
const word& name//gradf函数的第二个参数:场的名称,这个参数无关紧要
)
{
typedef typename outerProduct<vector, Type>::type GradType;//使用外积推导梯度场的类型,如果对一个标量场求梯度,那么会得到矢量场。如果对矢量场求梯度,那么会得到张量场,以此类推...
const fvMesh& mesh = ssf.mesh();//获取网格
tmp<GeometricField<GradType, fvPatchField, volMesh>> tgGrad//新建一个tmp临时变量,这个变量就是最终返回的梯度场。该场的数值是存储在网格体心的,而不是面心或者顶点
(
GeometricField<GradType, fvPatchField, volMesh>::New
(//GeometricField的构造函数,需要指定名称、网格、量纲、初始值、边界处数值外推方式
name,
mesh,
dimensioned<GradType>
(
"0",
ssf.dimensions()/dimLength,
Zero
),
extrapolatedCalculatedFvPatchField<GradType>::typeName
)
);
GeometricField<GradType, fvPatchField, volMesh>& gGrad = tgGrad.ref();//将tgGrad(tmp类型)引用赋值给gGrad(GeometricField)。之所以这样做是因为tmp类型仅仅用于储存对象而不能调用对象中的方法(函数),需要把对象的引用拿出来后才能调用其中的方法(函数)
const labelUList& owner = mesh.owner();//owner列表:内部面的所有者单元序号,来自const/polyMesh/owner文件的一部分
const labelUList& neighbour = mesh.neighbour();//neighbour列表:内部面的邻居单元序号,来自const/polyMesh/neighbour文件
const vectorField& Sf = mesh.Sf();//一个矢量场(场其实就是列表或数组),存储了内部面的法向量(不包含边界面),这些法向量的数值是存储在网格面心的,而不是体心或者顶点。法向量的大小是面的面积,方向由owner指向neighbour
Field<GradType>& igGrad = gGrad;//新建一个场,来自gGrad的内部网格(非边界),舍弃了gGrad的边界条件、量纲等信息。该场的数值是存储在网格体心的,而不是面心或者顶点。igGrad中的i表示internal,是内部的含义。
const Field<Type>& issf = ssf;//新建一个场,来自ssf的内部网格。该场的数值是存储在网格面心的,而不是体心或者顶点
//内部面对梯度的贡献
forAll(owner, facei)//遍历所有内部面(不包括边界面),面的序号是facei
{
GradType Sfssf = Sf[facei]*issf[facei];//计算面积与待求梯度场之乘积,Green-Gauss法
igGrad[owner[facei]] += Sfssf;//给第facei个面的owner单元增加Sfssf,因为Sf的方向由owner指向neighbour,朝向owner单元外部
igGrad[neighbour[facei]] -= Sfssf;//给第facei个面的neighbour单元减少Sfssf,因为Sf的方向由owner指向neighbour,朝向neighbour单元内部
}
//边界面对梯度的贡献
forAll(mesh.boundary(), patchi)//对每个patch遍历,注:往往边界是由若干个patch组成的
{
const labelUList& pFaceCells =
mesh.boundary()[patchi].faceCells();//一个列表,表示第patchi个patch所对应的单元序号,也就是这个patch上所有面所对应的owner单元序号。
const vectorField& pSf = mesh.Sf().boundaryField()[patchi];//一个向量场,表示第patchi个patch上的面法向量场。很显然,该场的数值是存储在网格面心的,而不是体心或者顶点。
const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];//边界上的ssf场
forAll(mesh.boundary()[patchi], facei)//对第patchi个patch,遍历每一个面
{
igGrad[pFaceCells[facei]] += pSf[facei]*pssf[facei];//给第facei个面的owner单元增加pSf[facei]*pssf[facei]
}
}
igGrad /= mesh.V();//igGrad除以单元体积即可得到梯度
gGrad.correctBoundaryConditions();//更新边界条件
return tgGrad;//返回tmp变量
}
openFoam中计算梯度都需要调用gradf函数,以fvc::grad函数为例,其返回值调用了gradf函数:
template<class Type>
tmp
<
GeometricField
<
typename outerProduct<vector, Type>::type, fvPatchField, volMesh
>
>
grad
(
const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
)
{
return fv::gaussGrad<Type>::gradf(ssf, "grad(" + ssf.name() + ')');
}
当然还有一系列grad函数的重载,但它们最终都回归到gradf函数