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{array}{}\text{(1)}& \frac{\mathrm{\partial }\rho }{\mathrm{\partial }t}=\mathrm{\nabla }\cdot \varphi \end{array}$$
明显二者量纲是不一致的呀。

[程迪]

@程迪

莫非OF中的div和数学$\mathrm{\nabla }\cdot$不一样，是这个意思：

$$\begin{array}{}\text{(2)}& \text{fvc::div}({\varphi }_f,{\alpha }_f)==\frac{\sum _f{\varphi }_f{\alpha }_f}{V}\end{array}$$

[程迪]

@程迪

还有一种可能是非定常项是乘以了体积V的。

[CFD中文网]

@程迪

$$\begin{array}{}\text{(3)}& \int \int \frac{\mathrm{\partial }\mathbf{U}}{\mathrm{\partial }t}\mathrm{d}V\mathrm{d}t=({\mathbf{U}}_{\mathrm{P}}^{n+1}-{\mathbf{U}}_{\mathrm{P}})V_{\mathrm{P}}\end{array}$$

没有细看你一楼的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{array}{}\text{(4)}& \text{fvc::ddt}(\alpha )=\frac{1}{\mathrm{\Delta }t}({\alpha }^n-{\alpha }^{n-1})\end{array}$$
这样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中文网]

嗯，没有显性处理时间项，而是在其他项除掉体积。

[CFD中文网]

@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]

[wwzhao]

@程迪 fvm下面的operator返回的是fvMatrix，而fvc下面的operator返回的是GeometricField。

fvMatrix是将物理量对单元体积分后离散得到的矩阵，由于对体积积分，所以需要将量纲乘以体积。

而之所以fvMatrix和GeometricField两种不同的类能相加相减等操作，是由于操作符已被重载。

[程迪]

@wwzhao

不太理解，你的意思是对于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是除以了体积的对角系数，这样$A{V}_{c}U$的量纲就能和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.C

tmp<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]

[wwzhao]

@程迪

对于fvm算子生成的fvMatrix和fvc生成的volScalarField相加生成fvMatrix的时候，volScalarField会被自动乘以网格体积？

没错，参考 https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C#L1691

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是除以了体积的对角系数

没错，参考 https://github.com/OpenFOAM/OpenFOAM-dev/blob/master/src/finiteVolume/fvMatrices/fvMatrix/fvMatrix.C#L724

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;
}