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两种不同的类能相加相减等操作,是由于操作符已被重载。



  • @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是除以了体积的对角系数,这样$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.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]


  • 版主

    @程迪

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

登录后回复
 

与 CFD中文网 的连接断开,我们正在尝试重连,请耐心等待