CFD Online Logo CFD Online URL
www.cfd-online.com
[Sponsors]
Home > Forums > OpenFOAM Programming & Development

Help!! customize surface tension term in interFoam

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

Reply
 
LinkBack Thread Tools Display Modes
Old   February 11, 2016, 16:59
Default Help!! customize surface tension term in interFoam
  #1
New Member
 
Join Date: May 2014
Posts: 5
Rep Power: 4
w051cxw is on a distinguished road
Hi guys, I am trying to replace the original surface tension term in 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 as a scalar and then reconstruct to the cell center as a vector. 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. 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.

Code:
    volVectorField newSurfaceTension //create volume vector field
    (
        IOobject
        (
            "newSurfaceTension",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::NO_WRITE
        ),
		mesh,
		dimensionedVector("newSurfaceTension",dimensionSet(1,-4,-2,0,0,0,0),Foam::vector(0,0,0))
    );

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
	{
		s=constant1;
	}
	else // if one is below 0.5 and the other is above 0.5
	{
		s=constant2;
	}
					
	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. Using UEqn.H as an example.
Code:
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      + turbulence->divDevRhoReff(rho, U)
    );

    UEqn.relax();

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                    //fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
                    newSurfaceTension2
                  - 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

Old   February 12, 2016, 03:06
Default
  #2
Senior Member
 
akidess's Avatar
 
Anton Kidess
Join Date: May 2009
Location: Delft, Netherlands
Posts: 1,163
Rep Power: 20
akidess will become famous soon enough
My guess is you calculate the new surface tension in createFields.H - it is never updated.
__________________
*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
akidess is offline   Reply With Quote

Old   February 12, 2016, 03:35
Default
  #3
New Member
 
Join Date: May 2014
Posts: 5
Rep Power: 4
w051cxw is on a distinguished road
Quote:
Originally Posted by akidess View Post
My guess is you calculate the new surface tension in createFields.H - it is never updated.

Thanks for your reply, Anton. I will try to implement the new surface tension to interFoam.C rather than createFields.H then. I will keep you updated.
w051cxw is offline   Reply With Quote

Old   February 12, 2016, 09:52
Default
  #4
Senior Member
 
chegdan's Avatar
 
Daniel P. Combest
Join Date: Mar 2009
Location: St. Louis, USA
Posts: 585
Rep Power: 20
chegdan will become famous soon enoughchegdan will become famous soon enough
Just to add more to this, the nested for loops on the order of n^2 may be slower than if you iterated over cell faces (instead of cell centers) and looked at owner and neighbour cells of that face. Look in the

$FOAM_SRC/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
__________________
Dan

Find me on twitter @dancombest and LinkedIn
chegdan is offline   Reply With Quote

Old   February 12, 2016, 15:55
Default
  #5
New Member
 
Join Date: May 2014
Posts: 5
Rep Power: 4
w051cxw is on a distinguished road
Hi Anton,

I tried to insert my piece of code (let's call it SURFACE_TENSION and can be found in my original post) to multiple places of interFoam.C (one at a time), however, the results are still like you said before, it was never updated.... The relevant code are below. Could anyone give me more advice?

Code:
    while (runTime.run())
    {
        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "alphaCourantNo.H"
        #include "setDeltaT.H"

        runTime++;

        Info<< "Time = " << runTime.timeName() << nl << endl;

        twoPhaseProperties.correct();

        #include "alphaEqnSubCycle.H"
        interface.correct();

        SURFACE_TENSION // first test is to insert here  
       
        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {

         SURFACE_TENSION // second test is to insert here

            #include "UEqn.H"

            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
            }
        }

        runTime.write();
w051cxw is offline   Reply With Quote

Old   February 12, 2016, 15:56
Default
  #6
New Member
 
Join Date: May 2014
Posts: 5
Rep Power: 4
w051cxw is on a distinguished road
Quote:
Originally Posted by chegdan View Post
Just to add more to this, the nested for loops on the order of n^2 may be slower than if you iterated over cell faces (instead of cell centers) and looked at owner and neighbour cells of that face. Look in the

$FOAM_SRC/finiteVolume/interpolation/surfaceInterpolation/limitedSchemes/LimitedScheme/LimitedScheme.C
Thanks Dan. I will go ahead and think about the calculation efficiency after my code works.
w051cxw is offline   Reply With Quote

Reply

Thread Tools
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 On
Pingbacks are On
Refbacks are On


Similar Threads
Thread Thread Starter Forum Replies Last Post
VOF with surface tension and wall adhesion Agad15 Fluent Multiphase 3 July 15, 2016 16:39
[ICEM] Problems with coedge curves and surfaces tommymoose ANSYS Meshing & Geometry 0 August 5, 2011 16:02
surface tension gradient on a free surface Abrem FLUENT 1 April 30, 2006 03:41
Surface tension - Continuum Surface Force model Robert Main CFD Forum 0 May 2, 2002 07:34
CFX4.3 -build analysis form Chie Min CFX 5 July 12, 2001 23:19


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