CFD Online Discussion Forums

CFD Online Discussion Forums (
-   OpenFOAM Running, Solving & CFD (
-   -   curvature at the interface (interFOAM) (

lindstroem August 27, 2012 05:36

curvature at the interface (interFOAM)
Hi all,

my question ist about the curvature calculation in interFoam. The magic happens in interfaceProperties.C and was also discussed often in the Forum. Anyhow, I still have an open question.
The curvature is calculated by K_ = -fvc::div(nHatf_); in the simple expression which is ok for me. Before that, the gradient of alpha is determined on the cell faces by interpolation from the cell-gradient of alpha. Afterwards, the face-unit interface normal flux is calculated (nHatf_ = nHatfv & Sf;) and used for the calculation of the curvature as written above.
But... I couldn't find out, how this divergence is realized. The Doxygen linked me to an explanation which is not helpful for me. I also set up a dummy hex 3x3 case an tried to understand the calculated value for the curvature in the central-cell, but did not succeed.

I'll be glad for any advices.

santiagomarquezd August 27, 2012 13:52

Hi Lindstroem, respect to the code of fvc::div, in gaussConvectionScheme.C we have:


00110 template<class Type>
00111 tmp<GeometricField<Type, fvPatchField, volMesh> >
00112 gaussConvectionScheme<Type>::fvcDiv
00113 (
00114    const surfaceScalarField& faceFlux,
00115    const GeometricField<Type, fvPatchField, volMesh>& vf
00116 ) const
00117 {
00118    tmp<GeometricField<Type, fvPatchField, volMesh> > tConvection
00119    (
00120        fvc::surfaceIntegrate(flux(faceFlux, vf))
00121    );

line 120 leads to fvcSurfaceIntegrate.C where after some jumps ends in:


00042 template<class Type>
00043 void surfaceIntegrate
00044 (
00045    Field<Type>& ivf,
00046    const GeometricField<Type, fvsPatchField, surfaceMesh>& ssf
00047 )
00048 {
00049    const fvMesh& mesh = ssf.mesh();
00051    const labelUList& owner = mesh.owner();
00052    const labelUList& neighbour = mesh.neighbour();
00054    const Field<Type>& issf = ssf;
00056    forAll(owner, facei)
00057    {
00058        ivf[owner[facei]] += issf[facei];
00059        ivf[neighbour[facei]] -= issf[facei];
00060    }
00062    forAll(mesh.boundary(), patchi)
00063    {
00064        const labelUList& pFaceCells =
00065            mesh.boundary()[patchi].faceCells();
00067        const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];
00069        forAll(mesh.boundary()[patchi], facei)
00070        {
00071            ivf[pFaceCells[facei]] += pssf[facei];
00072        }
00073    }
00075    ivf /= mesh.V();
00076 }

which basically implements (3.15) in H. PhD Thesis.


lindstroem August 28, 2012 05:37

Hi Santiago,

thanks for your help. But doesn't the GaussConvectionScheme need two parameters (faceFlux and vf)? in the interfaceProperties it is only called with one argument (nHatfv_).

The reason why I am digging through that, is that I wanted to calculate the sum of the curvature over a circle which should analytically be 1/r. If I sum up the K_ in the interfaceProperties it results in sth. about 10e-6 where the analytical solution would be 6.6. Do you know the reason for that?


akidess August 28, 2012 06:48

Right. You are calling fvc::div. You can see it also calls fvc::surfaceIntegrate.

lindstroem August 28, 2012 07:15

Hi Anton,

thanks, I think I got it :)


lindstroem August 29, 2012 09:12

Hi again,
I would like to specify the question posted above. I have a 2D case with a bubble. If I sum up the curvature in interfaceProperties.C

K_ = -fvc::div(nHatf_);

scalar KSum = 0.0;
forAll(gradAlpha, cells)
    KSum += K_[cells];
Info << KSum << endl;

I would expect to come "close" to the analytical solution of 1/r. But I am far away from that.
Has anyone did the same and can tell me what I did wrong?

akidess August 29, 2012 09:51

How did you initialize the shape?

lindstroem August 29, 2012 09:54

setFields with cylinderToCell

        p1 (0.5 0.5 -1);
        p2 (0.5 0.5 1);
        radius 0.25; 
            volScalarFieldValue alpha1 0
            volVectorFieldValue U (0 0 0)

The Domain is a square 1 x 1

akidess August 29, 2012 10:06

From what I know it's better to start with a non-equilibrium shape (e.g. a box), relax the shape and then evaluate the resulting curvature.

lindstroem August 29, 2012 10:18

Thanks for your comment. Just tried that:

The sum of the curvature starts with -6.802736e-11 and ends with -9.166001e-12. Should be something about 4.

//edit: seems to be known:

duongquaphim September 7, 2012 08:45

I don't understand why you want to check the accuracy of computed curvature in that way. Actually, at the interface of the droplet, at every point the curvature should be 1/r. I do not think summation of the interface curvature will equal 1/r. It should be n/r where n is the number of points in the summation.

Regarding to the curvature, you should note that the curvature vary quite a lot in interfoam. If you take the average of the curvature at the iso-surface of alpha1 = 0.5, you will get a value close to the value 1/r. Since the interface is smeared over 3-4 cells approximately, at cells with alpha1 = 0.9 or 0.1, you still have curvature different than 0. Therefore summing the curvature of cells belonging to the interface also can not give a proper result.

Hope it is clear for you.

lindstroem September 11, 2012 04:06

Hello Duong,

i got your point - absolutely right that summing up would result in sth like n/r. Thanks for your comment!


lindstroem September 19, 2012 10:46

Hi again,

I'd like to ask one more detail to my initial question concerning fvc::div():
Is it true, that we calculate the second derivative (the divergence) only at the centroid point of the cells using the first derivative (the gradients) from the faces?

akidess September 19, 2012 10:55

To me divergence and second derivative refer to two different things...

lindstroem September 19, 2012 11:07

ok, sorry, the divergence of the gradient...

All times are GMT -4. The time now is 00:12.