请教一个问题,关于Foam::fvMatrix<Type>::A() 函数的
-
template<class Type> Foam::tmp<Foam::volScalarField> Foam::fvMatrix<Type>::A() const { tmp<volScalarField> tAphi ( new volScalarField ( IOobject ( "A("+psi_.name()+')', psi_.instance(), psi_.mesh(), IOobject::NO_READ, IOobject::NO_WRITE ), psi_.mesh(), dimensions_/psi_.dimensions()/dimVol, zeroGradientFvPatchScalarField::typeName ) ); tAphi().internalField() = D()/psi_.mesh().V(); tAphi().correctBoundaryConditions(); return tAphi; }
为什么要除以
V()
???感觉跟李东岳的icoFoam解析对不上了。
tAphi().internalField() = D()/psi_.mesh().V();
-
试着回答一下。本身A()这个成员函数的目的就是返回代数矩阵的格心系数。而实际当中fvmatrix在构造时都已经包含了网格体积,即完成了体积分或者面积分运算。因此在取回格心系数(central coeffcient)时需要格外剔除网格体积项。否则在内网格部分就和D()这个成员函数重复了。
-
@youmengtian :kiss:
-
@youmengtian 你还记得在哪些文件或者代码里描述了这种处理的吗?如果可以的话,指点一下。谢谢
-
simpleFoam
-
在icoFoam中插入以下代码
fvVectorMatrix UEqn ( fvm::ddt(U) +fvm::div(phi,U) -fvm::laplacian(nu,U) ) Info<<"fvc::laplacian(nu,U).dimensions="<<fvc::laplacian(nu,U)().dimensions()<<endl; Info<<"fvm::laplacian(nu,U).dimensions="<<fvm::laplacian(nu,U)().dimensions()<<endl; Info<<"phi.dimensions="<<phi.dimensions()<<endl; Info<<"U.dimensions="<<U.dimensions()<<endl; Info<<"fvc::div(phi,U).dimensions="<<fvc::div(phi,U)().dimensions()<<endl; Info<<"fvm::div(phi,U).dimensions="<<fvm::div(phi,U)().dimensions()<<endl; Info << "Dimension of UEqn is "<<UEqn.dimensions()<<endl; // ...
输出大概是
fvc::laplacian(nu,U):[0 1 -2 0 0 0 0] fvm::laplacian(nu,U):[0 4 -2 0 0 0 0] phi:[0 3 -1 0 0 0 0] U:[0 1 -1 0 0 0 0] fvc::div(phi,U):[0 1 -2 0 0 0 0] fvm::div(phi,U):[0 4 -2 0 0 0 0] UEqn:[0 4 -2 0 0 0 0] UEqn.A():[0 0 -1 0 0 0 0] UEqn.H():[0 1 -2 0 0 0 0]
可以看出:
-
实际的情况是fvm项是乘以了体积的积分方程离散,输出类型是fvMatrix,fvc项是没有乘以体积的微分方程离散,输出类型是vol<Type>Field
-
从后面的代码
UEqn == -fvc::grad(p)
来看,二者是可以用operator==连接的。 -
根据fvMatrix的operator==源代码可以发现,operator==将GeometricField加到fvMatrix中的时候乘以了体积。
/* /src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C */ //line: 1458 template<class Type> Foam::tmp<Foam::fvMatrix<Type>> Foam::operator== ( const fvMatrix<Type>& A, const DimensionedField<Type, volMesh>& su ) { checkMethod(A, su, "=="); tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A)); tC.ref().source() += su.mesh().V()*su.field(); //!!! return tC; }
-
fvMatrix的部分函数比较诡异,UEqn.A()的量纲和UEqn相差[0 4 -1 0 0 0 0],而UEqn.A()和H()相差待求变量速度的量纲[0 1 -1 0 0 0 0],合理的解释是:
- UEqn是积分方程的离散,对于icoFoam,待求变量是速度U,量纲是[0 1 -1 0 0 0 0];
- UEqn方程系数和phi的量纲一致,是[0 3 -1 0 0 0 0],而方程右端源项是[0 4 -2 0 0 0 0]
- UEqn.A()是对角线系数除以体积,所以是[0 0 -1 0 0 0 0]
- UEqn.H()是非对角线部分作用于待求变量并被方程右端源项减去,再除以体积得到的,所以是[0 1 -2 0 0 0 0]
-
如果你详细看源代码,你会发现UEqn.D()是无量纲的,而UEqn.A()是强行把量纲附上去的,我觉得OF这些地方是有偷懒或者不一致的地方的。
-