CFD Online Logo CFD Online URL
Home > Forums > Software User Forums > OpenFOAM > OpenFOAM Running, Solving & CFD

Help!! customize surface tension term in interFoam

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

LinkBack Thread Tools Search this Thread Display Modes
Old   February 12, 2016, 02:15
Default Help!! customize surface tension term in interFoam
New Member
Join Date: May 2014
Posts: 5
Rep Power: 11
w051cxw is on a distinguished road
Hi guys, I am trying to replace the original surface tension term from interFoam using a custom code.

The original surface tension term mathematically can be expressed as Kappa*sigma*grad(alpha1). And interFoam expresses this term as fvc::interpolate(interface.sigmaK())*fvc::snGrad(a lpha1)

Based on my understanding, what interFoam does is to calculate the surface tension on face center and then reconstruct to the cell center. Therefore, "fvc::interpolate(interface.sigmaK())*fvc::snGrad( alpha1)" is a scalar rather than a vector.
What I am trying to do is to replace the above mentioned term by another function aiming to calculate the force between any cell to the rest of the cells in the domain.

Here is my code added to createFileds.H. It complies and runs OK, but no matter how long I run the code, the results of test case does not make sense. BTW, my test case is to place a air rectangle in a liquid box and see if the rectangle can become a circle eventually. Test case does not make sense means that the air rectangle does not change its shape.

volVectorField newSurfaceTension //create volume vector field

forAll(mesh.C(),celli) //calculate surface tension on celli

    forAll(mesh.C(),cellj) // relation between the celli to the rest of the cells in the domain

    if( ((alpha1[celli]<0.5)&& (alpha1[cellj]<0.5)) || ((alpha1[celli]>=0.5)&& (alpha1[cellj]>=0.5)) ) //if both less than 0.5 or both greater than 0.5
    else // if one is below 0.5 and the other is above 0.5
    newSurfaceTension[celli]=newSurfaceTension[celli]+(mesh.C()[celli]-mesh.C()[cellj])*s*constant3*mag(mesh.C()[celli]-mesh.C()[cellj]); // my formula


surfaceVectorField newSurfaceTension1("newSurfaceTension1",fvc::interpolate(newSurfaceTension)); //convert from volVectorField to surfaceVectorField
const surfaceVectorField& sf=mesh.Sf();
surfaceScalarField newSurfaceTension2("newSurfaceTension2",newSurfaceTension1 & sf); //covert from surfaceVectorField to surfaceScalarField
constant1, constant2, and constant3 are nothing but some constants.

After that I will update surface tension term in UEqn.H and pEqn.H. Here I use UEqn.H as an example.
    fvVectorMatrix UEqn
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      + turbulence->divDevRhoReff(rho, U)


    if (pimple.momentumPredictor())
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
I would greatly appreciate if anyone can give me any hint on why my code does not work.
w051cxw is offline   Reply With Quote


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
Access the surface tension force at a patch michielm OpenFOAM Programming & Development 1 December 13, 2013 04:50
Surface Tension term in two phase flows Soghrat Main CFD Forum 2 July 7, 2009 11:20
Surface tension asder CFX 2 August 4, 2008 10:41
Surface Tension term Soghrat Main CFD Forum 1 April 24, 2007 07:45
Surface tension - Continuum Surface Force model Robert Main CFD Forum 0 May 2, 2002 08:34

All times are GMT -4. The time now is 04:33.