[a.weber]

Hello everyone,

since i am working on a euler-lagrangian simulation of bubblecolumns, i stumbled over some papers concerning bubble (or particle) induce turbulence models.
One of them (seems to be the best fit to my physical system) is:
Pfleger, D., and S. Becker. "Modelling and simulation of the dynamic flow behaviour in a bubble column." Chemical Engineering Science 56.4 (2001): 1737-1747.

The equation to add is the following:


The idea is:
1. calculate the work/energy done by each particle on the fluid
"The energy input of the bubble wakes results from the forces acting between a gas bubble and the surrounding liquid and the local slip velocity." as shown in the paper.
Therefore I added a field "Work" similar to the UTrans calculation

In Foam::KinematicParcel<ParcelType>::calcVelocity:
....
// note: Feff.Sp() and Fc.Sp() must be the same
dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su());

// energy / work transported is needed for correct RANS simulation
dWork += (Ures.average() - Uc_) & dUTrans;

and

In Foam::KinematicParcel<ParcelType>::calc:
....
// Update momentum transfer coefficient
td.cloud().UCoeff()[cellI] += np0*Spu;

td.cloud().Work()[cellI] += np0*dWork;
....

2. Get this field and put it into the new kEpsilon Model by adding it to the "G" sourceterm:

volScalarField G(GName(),
    mut_*(tgradU() && dev(twoSymm(tgradU())))
    + C2_ * tWork //Line added from Pfleger Becker Equation
    );

Well, thats my overall idea.... still many problems.

Has anyone done this in OpenFOAM before or has ideas, how this could be done properly?
Main problem at the moment is to get the "Work" field into my turbulence part...

Greetings
Andy

[a.weber]

Update: Since the paper suggests to only take the drag-force into account, i will change to another solution:

I'm adding a drag vector to the particles variable list, which will allow me to also post-proccess the dragforce.
This list will then be used, to calculate my "work/energy" done to the fluid:

template<class CloudType>
inline const Foam::tmp<Foam::volScalarField>
Foam::KinematicCloud<CloudType>::work() const
{
    tmp<volScalarField> twork
    (
        new volScalarField
        (
            IOobject
            (
                this->name() + ":work",
                this->db().time().timeName(),
                this->db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            mesh_,
            dimensionedScalar("zero", dimPower, 0.0),
            zeroGradientFvPatchScalarField::typeName
        )
    );

    volScalarField& work = twork();
    forAllConstIter(typename KinematicCloud<CloudType>, *this, iter)
    {
        const parcelType& p = iter();
        const label cellI = p.cell();

        work[cellI] += p.nParticle()*mag(p.U()-p.Uc() & p.drag());
    }

    return twork;
}

hopefully, this will work soon.

[a.weber]

ready!
Now I have it working. A Pfleger Becker K Epsilon Turbulence Model.

If anyone is interested, just leave a reply.

[a.weber]

Hello again.
After some Testing, i found a little problem, when solving with mpirun.
The code runs without problems on a single node, but when it comes to multiprocessor runs I get an error.

As soon as the parcels are to be injected, my simulation will crash...only giving some default info:

[weber-VirtualBox:31670] *** An error occurred in MPI_Recv
[weber-VirtualBox:31670] *** on communicator MPI_COMM_WORLD
[weber-VirtualBox:31670] *** MPI_ERR_TRUNCATE: message truncated
[weber-VirtualBox:31670] *** MPI_ERRORS_ARE_FATAL: your MPI job will now abort

After some additional testing i found the problem sitting here:

template<class CloudType>
void Foam::InjectionModel<CloudType>::postInjectCheck
(
    const label parcelsAdded,
    const scalar massAdded
)
{
    const label allParcelsAdded = returnReduce(parcelsAdded, sumOp<label>()); // <--- bug with mpirun

    if (allParcelsAdded > 0)
    {
        Info<< nl
            << "--> Cloud: " << this->owner().name()
.....
    }

This is really confusing, because I didn't add anything to the code at this position. Still, it will always crash at this line. But why?

Does anyone know, what to do?

[Polli]

Hello Andi,
i have exactly the same problem like you with the injection problem in parallel.
Do you found any solution for the problem of the following error?

[weber-VirtualBox:31670] *** An error occurred in MPI_Recv [weber-VirtualBox:31670] *** on communicator MPI_COMM_WORLD [weber-VirtualBox:31670] *** MPI_ERR_TRUNCATE: message truncated [weber-VirtualBox:31670] *** MPI_ERRORS_ARE_FATAL: your MPI job will now abort

[a.weber]

Hi,

jep, this error was resolved in http://bugs.openfoam.org/view.php?id=1556

Just look at the OpenFOAM-2.4/dev version, there everything is working:
https://github.com/OpenFOAM/OpenFOAM...PatchInjection



Cheers
Andy

[suryakaundinya]

Hi,

I am currently trying to do the same thing and I have few questions. Did you add a new Parcel type or just made changes in Kinematic parcel directly and recompiled it? Because if you create a new parcel type the whole structure for the other parcel types has to be changed right?

With regards,
Surya

[a.weber]

I made a copy of src_intermediate and added the code to the kinematicParcel. Also made a copy of turbulenceModels and added stuff there.
Then just recompile with the modded src files.

[suryakaundinya]

Hello Andreas,

Can you please shed some light on how you called the "Work" field from the lagrangian model into the turbulence model.

Thanks and regards,
Surya.

[andyyy123]

Hi,

well, first add the function "work" like in my 1st post.
Then just initialize that as a field in the createClouds (or createFields):

work = parcels.work();
work.correctBoundaryConditions();

update this field after the lagrangian steps (in your main loop):

parcels.evolve();
...
work = parcels.work();
work.correctBoundaryConditions();

That should enable to read it in your turbulence model (or simply from whereever you have mesh access):

const tmp<volScalarField> work_ = mesh_.lookupObject<volScalarField>("work");

This is all similar to "alpha" or "theta" found in KinematicCloudI.H

Best regards
Andreas

[suryakaundinya]

Hi,
Thanks for the help. I did implement model as decribed above. But I have a small question (hope this will be the last one).
Did you define an initial condition for the source term seperately in 0/ folder during the simulation???
Surya.

[a.weber]

its set to uniform 0 at the start of every time step since you only want to have the work from the current step in your source.

[suryakaundinya]

Looks like I did something wrong. The simulation requires me to give an initial condition for this term.

[a.weber]

if I call parcel.work(), a temp field is created with values uniform 0.

template<class CloudType>
inline const Foam::tmp<Foam::volScalarField>
Foam::KinematicCloud<CloudType>::work() const
{
    tmp<volScalarField> twork
    (
        new volScalarField
        (
            IOobject
            (
                this->name() + ":work",
                this->db().time().timeName(),
                this->db(),
                IOobject::NO_READ,
                IOobject::NO_WRITE,
                false
            ),
            mesh_,
            dimensionedScalar("zero", dimPower/dimVol, 0.0),
            zeroGradientFvPatchScalarField::typeName
        )
);

also in my createField is use an empty field with NO_READ option

    volScalarField work
    (
        IOobject
        (
            "work",
            runTime.timeName(),
            mesh,
            IOobject::NO_READ,
            IOobject::AUTO_WRITE
        ),
        mesh,
        dimensionedScalar("zero", dimPower/dimVol, 0.0)
);

that should not ask for an initial condition

[suryakaundinya]

I did the same thing, but in the createFields.H I forgot to include
dimensionedScalar("zero", dimPower/dimVol, 0.0)
It's working now.
And one last clarification : Did you calculate the drag vector term directly in KinematicParcel.C from the SuSp and velocity terms?

[a.weber]

hmm not really....
unfortunately all coupled forces are summed up there.
so i save only the drag in an additional parcel property... (you could just take the "user" property for testing purpose)

i'm sure there is a better solution.... but this is what came to my mind at first

e.g. in breakup calculation
(void Foam::SprayParcel<ParcelType>::calcBreakup)
there may be a snippet of code you can use since there the "tmom" is calculated with forces:

    
const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, muAv);
const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, muAv);
this->tMom() = mass/(Fcp.Sp() + Fncp.Sp());

downside is!!! all forces are recalculated in this step (i don't know if it is better in v3.x or v4.x, i'm using v2.4)

[suryakaundinya]

Hi,

I did something similar. I created a function similar to calcVelocity, where I recalculated the forces.Later I used this function to calculate the drag force vector.

this->dragForce_= calcDrag(td, dt, cellI, Re, muc_, mass0);

Thanks for the inputs.
Surya.