# 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: 6 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: 418 Rep Power: 14 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. Post-doctoral Fellow Research Center for Computational Mechanics (CIMEC) - CONICET/FICH-UNL T.E.: 54-342-4511594 Ext. 1005 Güemes 3450 - (3000) Santa Fe Santa Fe - Argentina http://www.cimec.org.ar

 August 28, 2012, 05:37 #3 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 6 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: Delft, Netherlands Posts: 919 Rep Power: 17 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. *Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology

 August 28, 2012, 07:15 #5 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 6 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: 6 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: Delft, Netherlands Posts: 919 Rep Power: 17 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. *Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology

 August 29, 2012, 09:54 #8 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 6 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: Delft, Netherlands Posts: 919 Rep Power: 17 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. *Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology

 August 29, 2012, 10:18 #10 Senior Member   Join Date: Nov 2010 Posts: 113 Rep Power: 6 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: 92 Rep Power: 8 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: 6 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: 6 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: Delft, Netherlands Posts: 919 Rep Power: 17 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. *Help define the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oam-technology

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

 Thread Tools Display Modes Linear Mode

 Posting Rules You may not post new threads You may not post replies You may not post attachments You may not edit your posts BB code is On Smilies are On [IMG] code is On HTML code is OffTrackbacks are On Pingbacks are On Refbacks are On Forum Rules

 Similar Threads Thread Thread Starter Forum Replies Last Post Saturn CFX 34 October 16, 2014 05:27 duongquaphim OpenFOAM Running, Solving & CFD 23 July 25, 2013 01:03 Andrea_85 OpenFOAM 1 March 2, 2011 05:19 Pankaj CFX 9 November 23, 2009 05:05 lego CFX 3 November 5, 2002 21:09

All times are GMT -4. The time now is 09:14.