[Mahmoud_aboukhedr]

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

[hk318i]

For OpenFOAM-2.2;

First declare f as surfaceScalarField in createField.H

    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

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

fvc::interpolate(interface.sigmaK())*fvc::snGrad(alpha1)

, you can use

mixture.surfaceTensionForce()

Best wishes,
Hassan

[Mahmoud_aboukhedr]

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

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

    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;

    // 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

[hk318i]

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;

    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;

        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

    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;

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]

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:

    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:

    nHat_
    (
        IOobject
        (
            "nHat",
            gamma_.time().timeName(),
            gamma_.mesh()
        ),
    gamma_.mesh(),
    dimensionedVector("nHat", dimless, vector::zero)
    ),

In the interFoam solver, below the line

    interface.correct();

have:

volVectorField nHat = interface.nHat();
    if(runTime.write())
    {
        nHat.write();
    }

And inside createFields.H have:

    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.

[Mahmoud_aboukhedr]

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

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;

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,,