[gridley2]

Hello all,


I'm trying to modify some of the heat conduction terms in interCondensatingEvaporatingFoam in OpenFOAM-plus. I'm encountering some dimensions errors, and think I've created an example below which illustrates my confusion best. I expect that taking the divergence of the gradient of a dimensionless variable would have the dimensions of 1/length^2, but this seems to not be the case.


Consider this code:

 Info << "Units of alpha:" << alpha1.dimensions() << endl;                               
 surfaceScalarField mygrad(fvc::snGrad(alpha1));                                         
 Info << "Units of grad(alpha):" << mygrad.dimensions() << endl;                         
 volScalarField mydiv(fvc::div(mygrad));                                                 
Info << "Units of div(grad(alpha)):" << mydiv.dimensions() << endl;

Alpha1 is dimensionless here, so I expect to see:

Units of alpha:[0 0 0 0 0 0 0]
Units of grad(alpha):[0 -1 0 0 0 0 0]
Units of div(grad(alpha)):[0 -2 0 0 0 0 0]

Instead I see:

Units of alpha:[0 0 0 0 0 0 0]
Units of grad(alpha):[0 -1 0 0 0 0 0]
 Units of div(grad(alpha)):[0 -4 0 0 0 0 0]

What's up with that?


Calculating gradients at cell centers seems to give the expected result, using the below code:

    Info << "Units of alpha:" << alpha1.dimensions() << endl;
    volVectorField mygrad(fvc::grad(alpha1));
    Info << "Units of grad(alpha):" << mygrad.dimensions() << endl;
    volScalarField mydiv(fvc::div(mygrad));
    Info << "Units of div(grad(alpha)):" << mydiv.dimensions() << endl;

[massive_turbulence]

Originally Posted by gridley2

Hello all,


I'm trying to modify some of the heat conduction terms in interCondensatingEvaporatingFoam in OpenFOAM-plus. I'm encountering some dimensions errors, and think I've created an example below which illustrates my confusion best. I expect that taking the divergence of the gradient of a dimensionless variable would have the dimensions of 1/length^2, but this seems to not be the case.


Consider this code:

 Info << "Units of alpha:" << alpha1.dimensions() << endl;                               
 surfaceScalarField mygrad(fvc::snGrad(alpha1));                                         
 Info << "Units of grad(alpha):" << mygrad.dimensions() << endl;                         
 volScalarField mydiv(fvc::div(mygrad));                                                 
Info << "Units of div(grad(alpha)):" << mydiv.dimensions() << endl;

Alpha1 is dimensionless here, so I expect to see:

Units of alpha:[0 0 0 0 0 0 0]
Units of grad(alpha):[0 -1 0 0 0 0 0]
Units of div(grad(alpha)):[0 -2 0 0 0 0 0]

Instead I see:

Units of alpha:[0 0 0 0 0 0 0]
Units of grad(alpha):[0 -1 0 0 0 0 0]
 Units of div(grad(alpha)):[0 -4 0 0 0 0 0]

What's up with that?


Calculating gradients at cell centers seems to give the expected result, using the below code:

    Info << "Units of alpha:" << alpha1.dimensions() << endl;
    volVectorField mygrad(fvc::grad(alpha1));
    Info << "Units of grad(alpha):" << mygrad.dimensions() << endl;
    volScalarField mydiv(fvc::div(mygrad));
    Info << "Units of div(grad(alpha)):" << mydiv.dimensions() << endl;

Try this https://physics.stackexchange.com/qu...an-of-gaussian

Look in the comments!

[Daniel_Khazaei]

Originally Posted by massive_turbulence

Try this https://physics.stackexchange.com/qu...an-of-gaussian

Look in the comments!

well that also happens to keep my mind busy
The link you have posted is talking about divergence of gradient of a function called "g" that is of L^-2 dimension...then the dimension of the Laplacian of "g" becomes L^-4.

However, the "alpha" function in the first post is dimension-less...

[massive_turbulence]

Originally Posted by Daniel_Khazaei

well that also happens to keep my mind busy
The link you have posted is talking about divergence of gradient of a function called "g" that is of L^-2 dimension...then the dimension of the Laplacian of "g" becomes L^-4.

However, the "alpha" function in the first post is dimension-less...

I thought my link was also about dimensionless functions like sin and exp?

I'd like to see the actual proof myself.

[Daniel_Khazaei]

Originally Posted by massive_turbulence

I thought my link was also about dimensionless functions like sin and exp?

I'd like to see the actual proof myself.

well we need to find how Laplacian is defined in OpenFOAM...

[gridley2]

I appreciate the replies! Unfortunately, I understand that transcendental functions are dimensionless. The question here pertains to the fact that each spatial derivative should have units of inverse length, so a Laplacian should have inverse area units. Unfortunately it seems that obtaining gradients via snGrad then taking divergence gives a differing result (not inverse area, as it should be).



How exactly does that stack exchange post pertain to that?

[Daniel_Khazaei]

Originally Posted by gridley2

I appreciate the replies! Unfortunately, I understand that transcendental functions are dimensionless. The question here pertains to the fact that each spatial derivative should have units of inverse length, so a Laplacian should have inverse area units. Unfortunately it seems that obtaining gradients via snGrad then taking divergence gives a differing result (not inverse area, as it should be).

How exactly does that stack exchange post pertain to that?

I think we can find the answer here:

Meaning of "fvc::div(phi)"

[gridley2]

No, that doesn't explain the units here unfortunately. Even if you interpret divergence as integral of divergence over a cell volume, i.e. the sum of the surfaceScalarField times face areas, you should still not get the units which openfoam returns in the example above.

[Daniel_Khazaei]

Originally Posted by gridley2

No, that doesn't explain the units here unfortunately. Even if you interpret divergence as integral of divergence over a cell volume, i.e. the sum of the surfaceScalarField times face areas, you should still not get the units which openfoam returns in the example above.

This link I'm sure would explain what is going on here:

Incompatible dimensions....

fvc::div() actually divide the final result by the volume [m^3]

that's why the unit of div(phi) is [0 0 -1 0 0 0 0]

in your case: fvc::div(grad(Alpha)) -->> [0 0 -1 0 0 0 0] / [0 0 -3 0 0 0 0]


Edit: that was a mistake, I did somehow mix everything:
in your case: fvc::div(grad(Alpha)) -->> [0 0 -1 0 0 0 0] / [0 0 3 0 0 0 0]

[gridley2]

Thanks Daniel, I see what you're saying and agree that the post you provided gives some good insight. Unfortunately, it still doesn't make sense though; this does not give reason for divergence(grad(alpha)) to not have units of 1/length**2.


You say:

in your case: fvc::div(grad(Alpha)) -->> [0 0 -1 0 0 0 0] / [0 0 -3 0 0 0 0]

As a result of dividing by volume. Volume is units of [0 0 3 0 0 0 0], not [0 0 -3 0 0 0 0]. Also, since the grad(alpha) is defined as grad(alpha)*faceArea, the units should be [0 0 1 0 0 0 0].
Taking what they said about the fvc::div operator, you sum the surfaceScalarField values then divide by volume, this is [0 0 1 0 0 0 0] / [0 0 3 0 0 0 0] which should be units of 1/area i.e. [0 0 -2 0 0 0 0]. This is what I expect, but OpenFOAM clearly gives a different answer, as evidenced by the original post.

[Daniel_Khazaei]

Originally Posted by gridley2

Thanks Daniel, I see what you're saying and agree that the post you provided gives some good insight. Unfortunately, it still doesn't make sense though; this does not give reason for divergence(grad(alpha)) to not have units of 1/length**2.
You say:
As a result of dividing by volume. Volume is units of [0 0 3 0 0 0 0], not [0 0 -3 0 0 0 0]. Also, since the grad(alpha) is defined as grad(alpha)*faceArea, the units should be [0 0 1 0 0 0 0].
Taking what they said about the fvc::div operator, you sum the surfaceScalarField values then divide by volume, this is [0 0 1 0 0 0 0] / [0 0 3 0 0 0 0] which should be units of 1/area i.e. [0 0 -2 0 0 0 0]. This is what I expect, but OpenFOAM clearly gives a different answer, as evidenced by the original post.

I'm sorry that was a mistake when I copied the first dimension and I only changed [0 0 -1 0 0 0 0] to [0 0 -3 0 0 0 0]...(missed the -)!

What has been said on that post by "The surfaceField are presumed to be equal to field*Aface" is about "phi" as an example which is calculated by fvc::interpolate(U) & mesh.sf(), The unit of phi here will be:

[0 1 -1 0 0 0 0]*[0 2 0 0 0 0 0] = [0 3 -1 0 0 0 0]

Finally the sum of the resulting surfaceScalarField phi divided by the volume will be:

[0 3 -1 0 0 0 0] / [0 3 0 0 0 0 0] = [0 0 -1 0 0 0 0]

regarding your case, you know that fvc::snGrad(alpha1) already is a surfaceScalarField called "mygrad" which will have a dimension of:

[0 -1 0 0 0 0 0]

It is said that "the div operator is simply the sum of surfaceField on each cell-surfaces finally divided by the mesh volume"! So you shouldn't consider "mygrad*faceArea" as "mygrad" is a surfaceScalarField itself! By that I mean "fvc::div(mygrad)" sums the given surfaceScalarField "mygrad" over the cell faces which doesn't change the unit and then divide it by volume, so here we have:

[0 -1 0 0 0 0 0] / [0 3 0 0 0 0 0] = [0 -4 0 0 0 0 0]

That's my interpretation from that post...

[gridley2]

Thanks very much for the explanation Daniel! This makes sense. So to obtain the true divergence, I have to multiply "mygrad" by face areas first? To be honest, I'm surprised that ::div doesn't automatically do this. I wonder why not. I cannot thank you enough for helping me figure this out... I was certainly scratching my head over this for a few hours this past week!

[Daniel_Khazaei]

Originally Posted by gridley2

Thanks very much for the explanation Daniel! This makes sense. So to obtain the true divergence, I have to multiply "mygrad" by face areas first? To be honest, I'm surprised that ::div doesn't automatically do this. I wonder why not. I cannot thank you enough for helping me figure this out... I was certainly scratching my head over this for a few hours this past week!

well, what I think about and I'm not sure is:
surfaceScalarFields are variables defined on the faces of each cell and it should be used for variables that actually have physical meaning when defined on the surface, e.g. mass flux(rho*U*A)

fvc::div(anything) will give a result as long as you follow its rules and the provided argument doesn't violates any mathematical or C++ methods defined for it!

Whether the given argument and its divergence actually have any physical meaning or not is a different story and of-course non of fvc::div() business.