CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

curvature at the interface (interFOAM)

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   August 27, 2012, 06:36
Default curvature at the interface (interFOAM)
  #1
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 15
lindstroem is on a distinguished road
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
lindstroem is offline   Reply With Quote

Old   August 27, 2012, 14:52
Default
  #2
Senior Member
 
santiagomarquezd's Avatar
 
Santiago Marquez Damian
Join Date: Aug 2009
Location: Santa Fe, Santa Fe, Argentina
Posts: 452
Rep Power: 23
santiagomarquezd will become famous soon enough
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.
__________________
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
santiagomarquezd is offline   Reply With Quote

Old   August 28, 2012, 06:37
Default
  #3
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 15
lindstroem is on a distinguished road
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
lindstroem is offline   Reply With Quote

Old   August 28, 2012, 07:48
Default
  #4
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 29
akidess will become famous soon enough
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.
akidess is offline   Reply With Quote

Old   August 28, 2012, 08:15
Default
  #5
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 15
lindstroem is on a distinguished road
Hi Anton,

thanks, I think I got it

Greetings
lindstroem is offline   Reply With Quote

Old   August 29, 2012, 10:12
Default
  #6
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 15
lindstroem is on a distinguished road
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
lindstroem is offline   Reply With Quote

Old   August 29, 2012, 10:51
Default
  #7
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 29
akidess will become famous soon enough
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.
akidess is offline   Reply With Quote

Old   August 29, 2012, 10:54
Default
  #8
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 15
lindstroem is on a distinguished road
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
lindstroem is offline   Reply With Quote

Old   August 29, 2012, 11:06
Default
  #9
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 29
akidess will become famous soon enough
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.
akidess is offline   Reply With Quote

Old   August 29, 2012, 11:18
Default
  #10
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 15
lindstroem is on a distinguished road
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 09:11.
lindstroem is offline   Reply With Quote

Old   September 7, 2012, 09:45
Default
  #11
Member
 
Duong A. Hoang
Join Date: Apr 2009
Location: Delft, Netherlands
Posts: 93
Rep Power: 16
duongquaphim is on a distinguished road
Send a message via Yahoo to 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.
duongquaphim is offline   Reply With Quote

Old   September 11, 2012, 05:06
Default
  #12
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 15
lindstroem is on a distinguished road
Hello Duong,

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

Lindstroem
lindstroem is offline   Reply With Quote

Old   September 19, 2012, 11:46
Default
  #13
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 15
lindstroem is on a distinguished road
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?
lindstroem is offline   Reply With Quote

Old   September 19, 2012, 11:55
Default
  #14
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Germany
Posts: 1,377
Rep Power: 29
akidess will become famous soon enough
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.
akidess is offline   Reply With Quote

Old   September 19, 2012, 12:07
Default
  #15
Senior Member
 
Join Date: Nov 2010
Posts: 113
Rep Power: 15
lindstroem is on a distinguished road
ok, sorry, the divergence of the gradient...
lindstroem is offline   Reply With Quote

Reply

Thread Tools Search this Thread
Search this Thread:

Advanced Search
Display Modes

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 Off
Trackbacks are Off
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
Wind turbine simulation Saturn CFX 58 July 3, 2020 02:13
strange curvature with interFoam (comparison with Brackbill work) duongquaphim OpenFOAM Running, Solving & CFD 23 July 25, 2013 02: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 inlet-outlet pair lego CFX 3 November 5, 2002 21:09


All times are GMT -4. The time now is 22:35.