[lindstroem]

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.
Greetings
Lindstroem

[santiagomarquezd]

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();
00050 
00051     const labelUList& owner = mesh.owner();
00052     const labelUList& neighbour = mesh.neighbour();
00053 
00054     const Field<Type>& issf = ssf;
00055 
00056     forAll(owner, facei)
00057     {
00058         ivf[owner[facei]] += issf[facei];
00059         ivf[neighbour[facei]] -= issf[facei];
00060     }
00061 
00062     forAll(mesh.boundary(), patchi)
00063     {
00064         const labelUList& pFaceCells =
00065             mesh.boundary()[patchi].faceCells();
00066 
00067         const fvsPatchField<Type>& pssf = ssf.boundaryField()[patchi];
00068 
00069         forAll(mesh.boundary()[patchi], facei)
00070         {
00071             ivf[pFaceCells[facei]] += pssf[facei];
00072         }
00073     }
00074 
00075     ivf /= mesh.V();
00076 }

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

Regards.

[lindstroem]

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?

Greetings

[akidess]

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

[lindstroem]

Hi Anton,

thanks, I think I got it

Greetings

[lindstroem]

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?
Thanks!
Lindstroem

[akidess]

How did you initialize the shape?

[lindstroem]

setFields with cylinderToCell

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

The Domain is a square 1 x 1

[akidess]

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]

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: http://www.openfoam.org/mantisbt/pri...php?bug_id=158

[duongquaphim]

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]

Hello Duong,

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

Lindstroem

[lindstroem]

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]

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

[lindstroem]

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