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

Particle induced turbulence model - kEpsilon with additional source

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

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   July 21, 2015, 06:40
Default Particle induced turbulence model - kEpsilon with additional source
  #1
Member
 
Andreas Weber
Join Date: Jun 2014
Posts: 37
Rep Power: 10
a.weber is on a distinguished road
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:
G_{\alpha,BIT} = C_{k/\varepsilon} \mid M_{\alpha,\beta}\mid \cdot \mid u_\beta - u_\alpha\mid

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
Code:
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
Code:
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:
Code:
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 is offline   Reply With Quote

Old   July 27, 2015, 09:24
Default
  #2
Member
 
Andreas Weber
Join Date: Jun 2014
Posts: 37
Rep Power: 10
a.weber is on a distinguished road
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:

Code:
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.

Last edited by a.weber; July 28, 2015 at 06:22.
a.weber is offline   Reply With Quote

Old   July 28, 2015, 11:03
Default
  #3
Member
 
Andreas Weber
Join Date: Jun 2014
Posts: 37
Rep Power: 10
a.weber is on a distinguished road
ready!
Now I have it working. A Pfleger Becker K Epsilon Turbulence Model.

If anyone is interested, just leave a reply.
a.weber is offline   Reply With Quote

Old   July 30, 2015, 11:11
Default
  #4
Member
 
Andreas Weber
Join Date: Jun 2014
Posts: 37
Rep Power: 10
a.weber is on a distinguished road
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:
Code:
[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:
Code:
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?
a.weber is offline   Reply With Quote

Old   October 5, 2016, 04:30
Default
  #5
New Member
 
Werner
Join Date: Apr 2014
Posts: 19
Rep Power: 10
Polli is on a distinguished road
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
Polli is offline   Reply With Quote

Old   October 5, 2016, 08:49
Default
  #6
Member
 
Andreas Weber
Join Date: Jun 2014
Posts: 37
Rep Power: 10
a.weber is on a distinguished road
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
a.weber is offline   Reply With Quote

Old   April 24, 2017, 11:17
Default
  #7
New Member
 
surya kaundinya
Join Date: Mar 2017
Posts: 15
Rep Power: 7
suryakaundinya is on a distinguished road
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
suryakaundinya is offline   Reply With Quote

Old   April 28, 2017, 04:46
Default
  #8
Member
 
Andreas Weber
Join Date: Jun 2014
Posts: 37
Rep Power: 10
a.weber is on a distinguished road
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.
a.weber is offline   Reply With Quote

Old   August 1, 2017, 05:46
Default
  #9
New Member
 
surya kaundinya
Join Date: Mar 2017
Posts: 15
Rep Power: 7
suryakaundinya is on a distinguished road
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.
suryakaundinya is offline   Reply With Quote

Old   August 1, 2017, 11:52
Default
  #10
New Member
 
Andreas Weber
Join Date: Oct 2013
Posts: 1
Rep Power: 0
andyyy123 is on a distinguished road
Hi,

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

Code:
work = parcels.work();
work.correctBoundaryConditions();
update this field after the lagrangian steps (in your main loop):
Code:
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):
Code:
const tmp<volScalarField> work_ = mesh_.lookupObject<volScalarField>("work");
This is all similar to "alpha" or "theta" found in KinematicCloudI.H

Best regards
Andreas
andyyy123 is offline   Reply With Quote

Old   August 8, 2017, 10:37
Default
  #11
New Member
 
surya kaundinya
Join Date: Mar 2017
Posts: 15
Rep Power: 7
suryakaundinya is on a distinguished road
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.
suryakaundinya is offline   Reply With Quote

Old   August 8, 2017, 10:39
Default
  #12
Member
 
Andreas Weber
Join Date: Jun 2014
Posts: 37
Rep Power: 10
a.weber is on a distinguished road
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.
a.weber is offline   Reply With Quote

Old   August 8, 2017, 10:53
Default
  #13
New Member
 
surya kaundinya
Join Date: Mar 2017
Posts: 15
Rep Power: 7
suryakaundinya is on a distinguished road
Looks like I did something wrong. The simulation requires me to give an initial condition for this term.
suryakaundinya is offline   Reply With Quote

Old   August 8, 2017, 10:59
Default
  #14
Member
 
Andreas Weber
Join Date: Jun 2014
Posts: 37
Rep Power: 10
a.weber is on a distinguished road
if I call parcel.work(), a temp field is created with values uniform 0.
Code:
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
Code:
    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
a.weber is offline   Reply With Quote

Old   August 8, 2017, 11:12
Default
  #15
New Member
 
surya kaundinya
Join Date: Mar 2017
Posts: 15
Rep Power: 7
suryakaundinya is on a distinguished road
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?
suryakaundinya is offline   Reply With Quote

Old   August 9, 2017, 03:46
Default
  #16
Member
 
Andreas Weber
Join Date: Jun 2014
Posts: 37
Rep Power: 10
a.weber is on a distinguished road
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:
Code:
    
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)
a.weber is offline   Reply With Quote

Old   August 9, 2017, 06:17
Default
  #17
New Member
 
surya kaundinya
Join Date: Mar 2017
Posts: 15
Rep Power: 7
suryakaundinya is on a distinguished road
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.
suryakaundinya is offline   Reply With Quote

Reply

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
[swak4Foam] swak4foam building problem GGerber OpenFOAM Community Contributions 54 April 24, 2015 17:02
[swak4Foam] Swak4FOAM 0.2.3 / OF2.2.x installation error FerdiFuchs OpenFOAM Community Contributions 27 April 16, 2014 16:14
[swak4Foam] Error bulding swak4Foam sfigato OpenFOAM Community Contributions 18 August 22, 2013 13:41
OpenFOAM on MinGW crosscompiler hosted on Linux allenzhao OpenFOAM Installation 127 January 30, 2009 20:08
DecomposePar links against liblamso0 with OpenMPI jens_klostermann OpenFOAM Bugs 11 June 28, 2007 18:51


All times are GMT -4. The time now is 10:07.