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

Plotting Surface tension force with interfoam

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

Like Tree6Likes
  • 1 Post By hk318i
  • 2 Post By Mahmoud_aboukhedr
  • 3 Post By hk318i

Reply
 
LinkBack Thread Tools Display Modes
Old   September 11, 2015, 18:18
Default Plotting Surface tension force with interfoam
  #1
New Member
 
mahmoud aboukhedr
Join Date: Feb 2014
Location: London
Posts: 21
Rep Power: 5
Mahmoud_aboukhedr is on a distinguished road
ear Foamrs,

I was looking in all the threads for a clear way to plot the surface tension forces, but I could not find any.

So i tried to change the way interfoam (2.2.0) calculate the forces in interFoam,

So I changed the interface.C file and added this line

surfaceScalarField f = fvc::interpolate(interface.sigmaK())*fvc::snGrad(a lpha1);

Then I changed the Uequ.H File to be
fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(rhoPhi, U)
+ turbulence->divDevRhoReff(rho, U)
);

UEqn.relax();

if (pimple.momentumPredictor())
{
solve
(
UEqn
==
fvc::reconstruct
(
(
f
- ghf*fvc::snGrad(rho)
- fvc::snGrad(p_rgh)
) * mesh.magSf()
)
);
}

by doing so, no the solver calculate the surface tension force separately, then added to the momentum equation,
I tried it and it compiles and works fine...

Now I need to plot for the (f) and the interface curvature (K) in the processioning.

So I added these lines to the create field file;

volVectorField f
(
IOobject
(
"f",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimensionedVector("f",dimMass/(dimLength*dimLength*dimTime*dimTime), vector(0,0,0))
);



volScalarField K
(
IOobject
(
"K",
runTime.timeName(),
mesh,
IOobject::NO_READ,
IOobject::AUTO_WRITE
),
mesh,
dimless/dimTime
);

It compile fine .... but the problem that ;; when i plot for f or K ... the field is empty ..

So anyone can help ???
I can c f and K in each time step while am postpossessing .. but zero values ..

Any help will do ..

thanks
Mahmoud
Mahmoud_aboukhedr is offline   Reply With Quote

Old   September 12, 2015, 12:46
Default
  #2
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 231
Rep Power: 10
hk318i is on a distinguished road
For OpenFOAM-2.2;

First declare f as surfaceScalarField in createField.H

Code:
    surfaceScalarField f
    (
        IOobject
        (
            "f",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
    );
Then before the runTime.write() calculate f
Code:
f = fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1);
Finally compile the new solver.

Because it is surfaceField, it cannot loaded automatically in paraFoam. There is a small trick to load it as mentioned HERE

For OpenFOAM-2.3 and above, there are few changes in the transportModels library which give a new function to calculate the surface tension force directly.

Therefore instead of
Code:
fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
, you can use
Code:
mixture.surfaceTensionForce()
Best wishes,
Hassan
kiddmax likes this.
hk318i is offline   Reply With Quote

Old   September 12, 2015, 19:28
Default Plotting Surface tension force with interfoam
  #3
New Member
 
mahmoud aboukhedr
Join Date: Feb 2014
Location: London
Posts: 21
Rep Power: 5
Mahmoud_aboukhedr is on a distinguished road
Dear Hassen,

Thanks for your reply, that works fine for me. I will just re-post the steps for future reference so Foamers can use it.

So to plot surface scale field one can add the line in red to the end of the interFoam.C file as following
Code:
f = fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1);

        runTime.write();

        Info<< "ExecutionTime = " << runTime.elapsedCpuTime() << " s"
            << "  ClockTime = " << runTime.elapsedClockTime() << " s"
            << nl << endl;
    }

    Info<< "End\n" << endl;

    return 0;
Then When should change The Uequ.H with just adding f to the equation
Code:
    fvVectorMatrix UEqn
    (
        fvm::ddt(rho, U)
      + fvm::div(rhoPhi, U)
      + turbulence->divDevRhoReff(rho, U)
    );

    UEqn.relax();

    if (pimple.momentumPredictor())
    {
        solve
        (
            UEqn
         ==
            fvc::reconstruct
            (
                (
                  f
                  
                  - ghf*fvc::snGrad(rho)
                  - fvc::snGrad(p_rgh)
                ) * mesh.magSf()
            )
        );
    }
Then last thing to change the createfield.h file with adding this small code after calculating the alpha1 as following;

Code:
    // Construct interface from alpha1 distribution
    interfaceProperties interface(alpha1, U, twoPhaseProperties);



    surfaceScalarField f
    (
        IOobject
        (
            "f",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)
    );
Then to used in the post processing one can do the following;
1-Run foamToVTK -surfaceFields
2-Run paraFoam.
3-Load the surface fields base file "VTK/surfaceFields/surfaceFields_..vtk", so you can see them with the respective time snapshot.
4-Then apply the "Glyphs" filter to this file and you should see the respective points in glyph form.


Still ... am looking for plotting the Vol vector field !!
So any one can help in doing so??

Mahmoud
kmou and jamestangx like this.
Mahmoud_aboukhedr is offline   Reply With Quote

Old   September 13, 2015, 06:02
Default
  #4
Senior Member
 
Hassan Kassem
Join Date: May 2010
Location: UK
Posts: 231
Rep Power: 10
hk318i is on a distinguished road
If f will be used in UEqn, it should be calculated just before Pressure-velocity PIMPLE corrector loop.

To calculate f as volVectorField, define a new field in createFields as follows;

Code:
    volVectorField fV
    (
        IOobject
        (
            "fV",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::reconstruct(f* mesh.magSf())
    );
And in the main solver after updating f values, add;
Code:
        fV = fvc::reconstruct(f* mesh.magSf());
Memory-wise, I don't recommend this approach if you need only fV for post-processing. It is better to use the original UEqn and calculate only fV as follows;

in createFields
Code:
    volVectorField fV
    (
        IOobject
        (
            "fV",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        fvc::reconstruct((fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1))* mesh.magSf())
    );
In the main solver add;

Code:
fV = fvc::reconstruct((fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1))* mesh.magSf());
I think, there is still room for improvement but anyway it is better than calculating and writing f which is not needed.

Best wishes,
Hassan
kmou, akionux and han0459 like this.
hk318i is offline   Reply With Quote

Old   December 1, 2016, 08:42
Default
  #5
Member
 
Camille Bilger
Join Date: Jul 2013
Posts: 39
Rep Power: 5
kmou is on a distinguished road
Hi,


This all worked wonderfully, thank you.
I am trying to do a similar exercise with a different volVectorField nHat_ (interface unit normal) from interfaceProperties.C:

Code:
    volVectorField nHat_ = gradGamma/(mag(gradGamma) + deltaN_);
    forAll(nHat_.boundaryField(), patchi)
    {
        nHat_.boundaryField()[patchi] = nHatfv.boundaryField()[patchi];
    }

    K_ = -fvc::div(nHatf_) + (nHat_ & fvc::grad(nHatfv) & nHat_);

and define it in the constructor as:

Code:
    nHat_
    (
        IOobject
        (
            "nHat",
            gamma_.time().timeName(),
            gamma_.mesh()
        ),
    gamma_.mesh(),
    dimensionedVector("nHat", dimless, vector::zero)
    ),
In the interFoam solver, below the line
Code:
    interface.correct();
have:
Code:
volVectorField nHat = interface.nHat();
    if(runTime.write())
    {
        nHat.write();
    }
And inside createFields.H have:
Code:
    interfaceProperties interface(gamma, U, twoPhaseProperties);
    volVectorField nHat
    (
    IOobject
    (
        "nHatf",
        runTime.timeName(),
        mesh,
        IOobject::NO_READ,
        IOobject::AUTO_WRITE
    ),
    interface.nHat()
    );

But nHat does not seem to be passed to the solver correctly. Values are being printed inside interfaceProperties.C
But from the solver I get an empty internalField uniform(0 0 0) and at boundaries: calculated (0 0 0)

Any tips? It's strange that this procedure works for the volScalarField K but not for this volVectorField.
kmou is offline   Reply With Quote

Old   December 1, 2016, 09:37
Default Regarding nHat in post processing
  #6
New Member
 
mahmoud aboukhedr
Join Date: Feb 2014
Location: London
Posts: 21
Rep Power: 5
Mahmoud_aboukhedr is on a distinguished road
Hello Camille,

Am not 100% sure I got what you want to do, though I can see the problem.
For me the first part of the problem is using nHat for postprocessing like using K in postprocessing as we did before, If am not wrong.

For starter you can try adding those 4 lines to the main solver just to get nHat for post-processing;

In the main interFoam.C add
Code:
const volVectorField gradGamma = fvc::grad(gamma);
const dimensionedScalar deltaN("deltaN", dimensionSet(0,-1,0,0,0,0,0), 1E-16);
nHatF = (gradGamma/(mag(gradGamma)+deltaN));
Then In createFields.H add;
Code:
volVectorField nHatF
        (
            IOobject
            (
                "nHatF",
                runTime.timeName(),
                mesh,
                IOobject::MUST_READ,
                IOobject::AUTO_WRITE
            ),
            mesh
         );

        nHat.write();
and compile...

When you run the case, I will suppose it will stop and ask for nHatF file in the 0 folders >>
So just before running, make a copy of the U file and rename it nHatF ... and run
This should work fine just for getting nHat for post processing.
Let me know if it worked because it may not

If you want to use nhat for different calculation in the interfaceprop. I think you need to clarify a bit what exactly you want to do..

Also we still have the code master (hk318i) to check If he can help

Mahmoud,,
Mahmoud_aboukhedr 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
Problem with Gmsh nishant_hull Open Source Meshers: Gmsh, Netgen, CGNS, ... 23 August 5, 2015 02:09
[ICEM] Problems with coedge curves and surfaces tommymoose ANSYS Meshing & Geometry 0 August 5, 2011 16:02
Surface Tension Test Cases sega OpenFOAM Running, Solving & CFD 2 March 8, 2011 13:19
surface tension force Ray Main CFD Forum 0 October 17, 2001 15:42
gravitational force for free surface flow Jongtae Kim Main CFD Forum 1 July 2, 2000 11:57


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