CFD Online Discussion Forums (http://www.cfd-online.com/Forums/)
-   OpenFOAM Running, Solving & CFD (http://www.cfd-online.com/Forums/openfoam-solving/)
-   -   curvature at the interface (interFOAM) (http://www.cfd-online.com/Forums/openfoam-solving/106361-curvature-interface-interfoam.html)

 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.

Greetings
Lindstroem

 santiagomarquezd August 27, 2012 13:52

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

Code:

```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:

Code:

```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 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?

Greetings

 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 :)

Greetings

 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
Code:

```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 August 29, 2012 09:51

How did you initialize the shape?

 lindstroem August 29, 2012 09:54

setFields with cylinderToCell
Code:

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

 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

 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 13:19.