
[Sponsors] 
August 27, 2012, 05:36 
curvature at the interface (interFOAM)

#1 
Senior Member
Join Date: Nov 2010
Posts: 113
Rep Power: 7 
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 cellgradient of alpha. Afterwards, the faceunit 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 centralcell, 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: 430
Rep Power: 16 
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 ); 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 } Regards.
__________________
Santiago MÁRQUEZ DAMIÁN, Ph.D. Research Scientist Research Center for Computational Methods (CIMEC)  CONICET/UNL Tel: 543424511594 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: 7 
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 10e6 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: 1,141
Rep Power: 20 
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. *Join the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oHPxcPqde7HtA2 

August 28, 2012, 07:15 

#5 
Senior Member
Join Date: Nov 2010
Posts: 113
Rep Power: 7 
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: 7 
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; 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: 1,141
Rep Power: 20 
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. *Join the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oHPxcPqde7HtA2 

August 29, 2012, 09:54 

#8 
Senior Member
Join Date: Nov 2010
Posts: 113
Rep Power: 7 
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) ); 

August 29, 2012, 10:06 

#9 
Senior Member
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 1,141
Rep Power: 20 
From what I know it's better to start with a nonequilibrium 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. *Join the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oHPxcPqde7HtA2 

August 29, 2012, 10:18 

#10 
Senior Member
Join Date: Nov 2010
Posts: 113
Rep Power: 7 
Thanks for your comment. Just tried that:
The sum of the curvature starts with 6.802736e11 and ends with 9.166001e12. 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

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 isosurface of alpha1 = 0.5, you will get a value close to the value 1/r. Since the interface is smeared over 34 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: 7 
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: 7 
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: 1,141
Rep Power: 20 
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. *Join the OpenFOAM stackexchange Q&A site: http://area51.stackexchange.com/prop...oHPxcPqde7HtA2 

September 19, 2012, 11:07 

#15 
Senior Member
Join Date: Nov 2010
Posts: 113
Rep Power: 7 
ok, sorry, the divergence of the gradient...


Thread Tools  
Display Modes  


Similar Threads  
Thread  Thread Starter  Forum  Replies  Last Post 
Wind turbine simulation  Saturn  CFX  45  February 8, 2016 05:42 
strange curvature with interFoam (comparison with Brackbill work)  duongquaphim  OpenFOAM Running, Solving & CFD  23  July 25, 2013 01:03 
Curvature of the interface using interFoam  Andrea_85  OpenFOAM  1  March 2, 2011 05:19 
RPM in Wind Turbine  Pankaj  CFX  9  November 23, 2009 05:05 
Replace periodic by inletoutlet pair  lego  CFX  3  November 5, 2002 21:09 