fvc::div的前后量纲关系到底是啥？
-
如题，参考rhoCentralFoam.C代码178-179行。
// --- Solve density solve(fvm::ddt(rho) + fvc::div(phi));phi变量根据前面的定义是
phi = aphiv_pos*rho_pos + aphiv_neg*rho_neg; //159行：[phi]= 密度*[aphiv_pos] surfaceScalarField aphiv_pos("aphiv_pos", phiv_pos - aSf);//136行：[aphiv_pos]= [phiv_pos ] surfaceScalarField phiv_pos("phiv_pos", U_pos & mesh.Sf());//93行：[phiv_pos]=速度*面积phi的量纲是M/T，也就是质量流量(mass flow rate, [M/T])，而不是单位面积的质量流率(mass flux, [M/L^2/T])
而密度方程的意思应该是
\begin{equation}
\frac{\partial \rho}{\partial t}=\nabla\cdot\phi
\end{equation}
明显二者量纲是不一致的呀。 -
莫非OF中的div和数学$\nabla\cdot$不一样，是这个意思：
\begin{equation}
\text{fvc::div}(\phi_f,\alpha_f) == \frac{\sum_f{\phi_f\alpha_f}}{V}
\end{equation} -
还有一种可能是非定常项是乘以了体积V的。
-
\begin{equation}
\int\int\frac{\partial\mathbf{U}}{\partial t}\mathrm{d}V\mathrm{d}t=(\mathbf{U}_\mathrm{P}^{n+1}-\mathbf{U} _\mathrm{P})V _\mathrm{P}
\end{equation}没有细看你一楼的
alphav什么的，不过从方程离散来看，你最后的解释是正确的。公式显示有问题？我处理一下，
我处理好了:sunglasses:
-
@cfd-china
我觉得你是错的，应该是
以Euler格式为例:/*src\finiteVolume\finiteVolume\ddtSchemes\EulerDdtScheme\EulerDdtScheme.C*/ //line: 102 template<class Type> tmp<GeometricField<Type, fvPatchField, volMesh>> EulerDdtScheme<Type>::fvcDdt ( const GeometricField<Type, fvPatchField, volMesh>& vf ) { dimensionedScalar rDeltaT = 1.0/mesh().time().deltaT(); IOobject ddtIOobject ( "ddt("+vf.name()+')', mesh().time().timeName(), mesh() ); if (mesh().moving()) { return tmp<GeometricField<Type, fvPatchField, volMesh>> ( new GeometricField<Type, fvPatchField, volMesh> ( ddtIOobject, rDeltaT* ( vf() - vf.oldTime()()*mesh().Vsc0()/mesh().Vsc() //动网格也是同理 ), rDeltaT.value()* ( vf.boundaryField() - vf.oldTime().boundaryField() ) ) ); } else { return tmp<GeometricField<Type, fvPatchField, volMesh>> ( new GeometricField<Type, fvPatchField, volMesh> ( ddtIOobject, rDeltaT*(vf - vf.oldTime()) //显然如果vf的量纲是D,那么ddt(vf)的量纲是D/T ) ); } }翻译成公式是：
\begin{equation}
\text{fvc::ddt}(\alpha)=\frac{1}{\Delta t}(\alpha^n-\alpha^{n-1})
\end{equation}
这样OpenFOAM的代码中的公式量纲才是正确的。 -
@cfd-china
确认了，OpenFOAM嘴上说的和手底下做的不一致。
/*src/finiteVolume/finiteVolume/fvc/fvcSurfaceIntegrate.H*/ //line:79 template<class Type> tmp<GeometricField<Type, fvPatchField, volMesh>> surfaceIntegrate ( const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf ) { const fvMesh& mesh = ssf.mesh(); tmp<GeometricField<Type, fvPatchField, volMesh>> tvf ( new GeometricField<Type, fvPatchField, volMesh> ( IOobject ( "surfaceIntegrate("+ssf.name()+')', ssf.instance(), mesh, IOobject::NO_READ, IOobject::NO_WRITE ), mesh, dimensioned<Type> ( "0", ssf.dimensions()/dimVol, //**HERE!!!!!!!!** Zero ), extrapolatedCalculatedFvPatchField<Type>::typeName ) ); GeometricField<Type, fvPatchField, volMesh>& vf = tvf.ref(); surfaceIntegrate(vf.primitiveFieldRef(), ssf); vf.correctBoundaryConditions(); return tvf; } -
嗯，没有显性处理时间项，而是在其他项除掉体积。
-
@cfd-china 同样的，时间项除掉了时间，也就是本贴公式中的公式(4):joking: 。
@李东岳 或许可以更新一下icoFoam解析，目前OpenFOAM这种颠倒的操作，可能会引起困惑。 -
@cfd-china
今儿又发现一个吊诡的地方：fvm::div(phi,X)的量纲是什么？似乎TNND和fvc::div又不一样...
/* OpenFOAM-dev/src/finiteVolume/finiteVolume/convectionSchemes/gaussConvectionScheme/gaussConvectionScheme.C */ //line: 90 tmp<fvMatrix<Type>> tfvm ( new fvMatrix<Type> ( vf, faceFlux.dimensions()*vf.dimensions() ) );有人能解释一下么？
-
@cfd-china
关于OF中icoFoam的量纲测试
// 在合适的地方插入以下代码 fvVectorMatrix UEqn ( fvm::ddt(U) +fvm::div(phi,U) -fvm::laplacian(nu,U) ) Info << "Dimension of UEqn is "<<UEqn.dimensions()<<endl;输出显示量纲为
Dimension of UEqn is [0 4 -2 0 0 0 0]
相比之下，速度U的量纲为[0 1 -1 0 0 0 0]，体积流量phi的量纲是[0 3 -1 0 0 0 0]
所以实际的情况似乎又是fvm::ddt项乘以了体积，fvm::laplacian项乘以了体积，fvm::div项没有乘也没有除以体积。可能我前面的代码找到的地方不是fvm::ddt最终实现的地方吧。
注：icoFoam中的压力p是kinematic pressure. 量纲是压力除以密度[0 2 -2 0 0 0 0]
-
@程迪 fvm下面的operator返回的是fvMatrix，而fvc下面的operator返回的是GeometricField。
fvMatrix是将物理量对单元体积分后离散得到的矩阵，由于对体积积分，所以需要将量纲乘以体积。
而之所以fvMatrix和GeometricField两种不同的类能相加相减等操作，是由于操作符已被重载。
-
不太理解，你的意思是对于fvm算子生成的fvMatrix和fvc生成的volScalarField相加生成fvMatrix的时候，volScalarField会被自动乘以网格体积？
我觉得也有这种可能，我对icoFoam里面fvMatrix::A()和H()也输出dimensions看了一下。
rAU=1.0/UEqn.A()的量纲是[0 0 1 0 0 0 0]，所以UEqn.A()的量纲是[0 0 -1 0 0 0 0]，暂，感觉A是除以了体积的对角系数，这样$AV_cU$的量纲就能和UEqn的量纲对上了。HbyA的量纲是[0 1 -1 0 0 0 0]，和速度一样。所以H的量纲似乎应该是[0 1 -2 0 0 0 0]，也是除以过体积的非对角项之和，然后phiHbyA的量纲是[0 3 -1 0 0 0 0]。
比较诡异的是H的量纲和UEqn的量纲是不一样的，OF这事儿办得忒不地道了。
-
@cfd-china 多谢，有空更新一下。
@程迪 OpenFOAM这方面操作很绕。你输出的是矩阵的单位。以时间项为例：
EulerDdtScheme.Ctmp<fvMatrix<Type> > tfvm ( new fvMatrix<Type> ( vf, vf.dimensions()*dimVol/dimTime//fvMatrix单位设定 ) );比如我们速度
fvm:ddt(U)，离散后矩阵单位后：速度单位/时间单位×体积单位为：Dimension of UEqn is [0 4 -2 0 0 0 0] -
对于fvm算子生成的fvMatrix和fvc生成的volScalarField相加生成fvMatrix的时候，volScalarField会被自动乘以网格体积？
template<class Type> Foam::tmp<Foam::fvMatrix<Type>> Foam::operator+ ( const fvMatrix<Type>& A, const tmp<GeometricField<Type, fvPatchField, volMesh>>& tsu ) { checkMethod(A, tsu(), "+"); tmp<fvMatrix<Type>> tC(new fvMatrix<Type>(A)); tC.ref().source() -= tsu().mesh().V()*tsu().primitiveField(); tsu.clear(); return tC; }感觉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, extrapolatedCalculatedFvPatchScalarField::typeName ) ); tAphi.ref().primitiveFieldRef() = D()/psi_.mesh().V(); tAphi.ref().correctBoundaryConditions(); return tAphi; }
