# curvature at the interface (interFOAM)

 Register Blogs Members List Search Today's Posts Mark Forums Read

 August 27, 2012, 05:36 curvature at the interface (interFOAM) #1 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 15 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

 August 27, 2012, 13:52 #2 Senior Member     Santiago Marquez Damian Join Date: Aug 2009 Location: Santa Fe, Santa Fe, Argentina Posts: 452 Rep Power: 23 Hi Lindstroem, respect to the code of fvc::div, in gaussConvectionScheme.C we have: Code: ```00110 template 00111 tmp > 00112 gaussConvectionScheme::fvcDiv 00113 ( 00114 const surfaceScalarField& faceFlux, 00115 const GeometricField& vf 00116 ) const 00117 { 00118 tmp > tConvection 00119 ( 00120 fvc::surfaceIntegrate(flux(faceFlux, vf)) 00121 );``` line 120 leads to fvcSurfaceIntegrate.C where after some jumps ends in: Code: ```00042 template 00043 void surfaceIntegrate 00044 ( 00045 Field& ivf, 00046 const GeometricField& 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& 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& 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. __________________ Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC) - CONICET/UNL Tel: 54-342-4511594 Int. 7032 Colectora Ruta Nac. 168 / Paraje El Pozo (3000) Santa Fe - Argentina. http://www.cimec.org.ar

 August 28, 2012, 05:37 #3 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 15 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

 August 28, 2012, 06:48 #4 Senior Member     Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,377 Rep Power: 30 Right. You are calling fvc::div. You can see it also calls fvc::surfaceIntegrate. __________________ *On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer.

 August 28, 2012, 07:15 #5 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 15 Hi Anton, thanks, I think I got it Greetings

 August 29, 2012, 09:12 #6 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 15 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

 August 29, 2012, 09:51 #7 Senior Member     Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,377 Rep Power: 30 How did you initialize the shape? __________________ *On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer.

 August 29, 2012, 09:54 #8 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 15 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

 August 29, 2012, 10:06 #9 Senior Member     Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,377 Rep Power: 30 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. __________________ *On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer.

 August 29, 2012, 10:18 #10 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 15 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 Last edited by lindstroem; September 7, 2012 at 08:11.

 September 7, 2012, 08:45 #11 Member   Duong A. Hoang Join Date: Apr 2009 Location: Delft, Netherlands Posts: 93 Rep Power: 17 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.

 September 11, 2012, 04:06 #12 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 15 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 #13 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 15 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?

 September 19, 2012, 10:55 #14 Senior Member     Anton Kidess Join Date: May 2009 Location: Germany Posts: 1,377 Rep Power: 30 To me divergence and second derivative refer to two different things... __________________ *On twitter @akidTwit *Spend as much time formulating your questions as you expect people to spend on their answer.

 September 19, 2012, 11:07 #15 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 15 ok, sorry, the divergence of the gradient...