[Mehdi3031]

Hello All, I realized something really weird in OpenFoam which I can't find any explanation...
consider the solver interFoam, in creatField.H phi is made as incompressible flow:
#include "createPhi.H"
which in the definition can be found that:

surfaceScalarField phi
(
    IOobject
    (
        "phi",
        runTime.timeName(),
        mesh,
        IOobject::READ_IF_PRESENT,
        IOobject::AUTO_WRITE
    ),
    linearInterpolate(U) & mesh.Sf()
);

Now, U is [m/s] and Sf is [m^2], so phi should be [m^3/s]. If you printout the dimension of phi from the solver, you'll get the same answer:

    Info << "phi.dim = "<< phi.dimensions() << endl;
result:
phi.dim = [0 3 -1 0 0 0 0]

SO FAR SO GOOD!!!!
now lets consider the UEqn.H:

    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
	+ fvm::div(rhoPhi, U)
        + turbulence->divDevRhoReff(U)
    );

which I rewrite it as following so that we can see it better:

    volTensorField SS=fvc::grad(U);
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
	+ fvm::div(rhoPhi, U)
	- fvm::laplacian(mu, U)
	- fvc::div(  mu*dev2(SS.T())   )
    );

1)Dimension analysis of term "fvm::ddt(rho, U)":
[1/s]*[kg/m^3]*[m/s]= [kg/m^2/s^2] ([1/s] is due to derivative with respect to time)

2)Dimension analysis of term "fvm::laplacian(mu, U)":
[1/m]*[1/m]*[kg/m/s]*[m/s]=[kg/m^2/s^2] ([1/m] is due to derivative with respect to space , "gradient" )

so the divergence term should also have the same unit, lets see:

3)Dimension analysis of term "fvm::div(rhoPhi, U)":
[1/m]*[kg/m^3]*[m^3/s]*[m/s]=[kg/s^2] ([1/m] is due to derivative with respect to space , "gradient" )


WHYYYYYYYYY????

now if u check the last term, you'll realize that the only way to make this term to have the same unit as other term, phi should be [m/s].

ANY COMMENT??

[Mehdi3031]

by the way, this is the case for any other solvers as well!! Not just interFoam...

[Zeppo]

Don't forget that fvMatrix equations are not formulated in a derivative form but rather in an integral form.It means that every term in your matrix equation represents an integral not derivative. Yes, initially an equation is formulated in a derivative form and then it is integrated over the control volume. Applying Gauss theorem some terms (div, laplacian) are converted from volume integral to surface integral. To learn more about the discretization of equations in FVM see books like these:

1. Moukalled, Mangani, Darwish - The Finite Volume Method in Computational Fluid Dynamics An Advanced Introduction with OpenFOAM® and Matlab (Fluid Mechanics and Its Applications) - 2015
2. Anderson J. - Computational fluid dynamics. The basics with applications (MGH, 1995)(T)(563s)
3. Versteeg H, Malalasekera W - An Introduction to Computational Fluid Dynamics. The Finite Volume Method - 2007
4. Ferziger J.H., Peric M. Computational Methods for Fluid Dynamics - 2002

Source code (as always) is your friend. I marked the relevant lines in code with red color.

1) fvm::ddt(rho, U):
[kg/m^3]*[m/s]*[m^3]/[s] = [kg*m/s^2]

template<class Type>
tmp<fvMatrix<Type> >
EulerDdtScheme<Type>::fvmDdt
(
    const volScalarField& rho,
    const GeometricField<Type, fvPatchField, volMesh>& vf
)
{
    tmp<fvMatrix<Type> > tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            rho.dimensions()*vf.dimensions()*dimVol/dimTime
        )
    );
    fvMatrix<Type>& fvm = tfvm();

    scalar rDeltaT = 1.0/mesh().time().deltaTValue();

    fvm.diag() = rDeltaT*rho.internalField()*mesh().Vsc();

    if (mesh().moving())
    {
        fvm.source() = rDeltaT
            *rho.oldTime().internalField()
            *vf.oldTime().internalField()*mesh().Vsc0();
    }
    else
    {
        fvm.source() = rDeltaT
            *rho.oldTime().internalField()
            *vf.oldTime().internalField()*mesh().Vsc();
    }

    return tfvm;
}

3) fvm::div(rhoPhi, U):
[kg/m^3]*[m^3/s]*[m/s]=[kg*m/s^2]

template<class Type>
tmp<fvMatrix<Type> >
gaussConvectionScheme<Type>::fvmDiv
(
    const surfaceScalarField& faceFlux,
    const GeometricField<Type, fvPatchField, volMesh>& vf
) const
{
    tmp<surfaceScalarField> tweights = tinterpScheme_().weights(vf);
    const surfaceScalarField& weights = tweights();

    tmp<fvMatrix<Type> > tfvm
    (
        new fvMatrix<Type>
        (
            vf,
            faceFlux.dimensions()*vf.dimensions()
        )
    );
    fvMatrix<Type>& fvm = tfvm();

    fvm.lower() = -weights.internalField()*faceFlux.internalField();
    fvm.upper() = fvm.lower() + faceFlux.internalField();
    fvm.negSumDiag();

    forAll(vf.boundaryField(), patchI)
    {
        const fvPatchField<Type>& psf = vf.boundaryField()[patchI];
        const fvsPatchScalarField& patchFlux = faceFlux.boundaryField()[patchI];
        const fvsPatchScalarField& pw = weights.boundaryField()[patchI];

        fvm.internalCoeffs()[patchI] = patchFlux*psf.valueInternalCoeffs(pw);
        fvm.boundaryCoeffs()[patchI] = -patchFlux*psf.valueBoundaryCoeffs(pw);
    }

    if (tinterpScheme_().corrected())
    {
        fvm += fvc::surfaceIntegrate(faceFlux*tinterpScheme_().correction(vf));
    }

    return tfvm;
}

[Mehdi3031]

Bravooo, I see now, Thank you so much for your clear guide