Two way particle coupling in OF 2.0 and the particle reynolds number
 

Hi All,

I have been trying to get two way particle coupling working in OF 2.0.
Essentially I have modified as required Bernhard Gschaider's icoLagrangianFoam (which works well in OF 1.6-ext) so that it will also (hopefully) work in OF 2.0.
I think that I more or less have it working (well the particles and the bulk flow are coupled), but the particles seem to slow down much faster than they should.

To try and track this down I have been mapping out various variables to see where things seem to have gone wrong. In my test case I am assuming that the particles are spheres, and so they should use the sphere drag model.
Anyway when I map out (print to the screen) the Reynolds number used in calculating the Drag Coefficient I am finding that the Reynolds number is being calculated using the particle density as opposed to the bulk fluid density (hence the Reynolds number is larger than it should be and the drag is then also significantly larger, which would account for why the particles slow down faster than they should).
I have cross referenced this with the same case running in OF 1.6-ext and I can confirm that it uses the correct Reynolds number.

As the largrangian source terms returned by OF 2.0 seem to be Forces (as opposed to acceleration terms), I have needed to multiply the equations through by rho (so this could possible explain where the error lies, although as the test problem has water particles in air it should only be out by the air density and not the water density (although I could be missing something here)).

So my questions are essentially two fold, firstly has anyone seen this before? and secondly where can I find (name of the routine or header file) where the Particle Reynolds number (for the purposes of calculating the particles drag coefficient) is calculated? (so I can print out the values and check that they are what I think they should be).

Thanks in advance,

Paul.

Hi,

The SourceTerm has the dimension of a force and u need to adapt the SourceTerm for the NSE.
Obviously, if u start the calculation just with the SourceTerm SU() u get an error, because the dimensions dont match.
Take a look at the code for SU() in the KinematicCloudI.H.

U can find the calculation for Re in KinematicParcelI.H. (http://foam.sourceforge.net/docs/cpp...ce.html#l00534)
The calculation for the DragCoeff (multiplied by the Re-number) is in SphereDragForce.C (submodels -> kinematic -> particleforces) -> CdRe(Re)

Regards,
Jens

Hi Jens,

Thanks for your reply. The source code I am using is listed below. I used the info from KinematicCloudI.H to modify the solver as I thought necessary, I also needed to modify a couple of the header files (only minor changes).

If anyone can spot any obvious errors please let me know.

Thanks, I will have a look at how Re is calculated in KinematicParcelI.H.. I was monitoring the Re number by writing out the value used in the CdRe calculation in SphereDragForce.C. So that was where I was observing that the Re values that were being reported seemed to be calculated using the particle density instead of the fluid density.

Thanks,

Paul.

\*---------------------------------------------------------------------------*/

#include "fvCFD.H"
#include "singlePhaseTransportModel.H"
#include "turbulenceModel.H"
#include "basicKinematicCollidingCloud.H"
#include "pimpleControl.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

int main(int argc, char *argv[])
{
argList::addOption
(
"cloudName",
"name",
"specify alternative cloud name. default is 'kinematicCloud'"
);

#include "setRootCase.H"
#include "createTime.H"
#include "createMesh.H"
#include "createFields.H"
#include "initContinuityErrs.H"

// * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * * //

Info<< "\nStarting time loop\n" << endl;

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

Info<< "Evolving " << kinematicCloud.name() << endl;

# include "readPISOControls.H"
# include "CourantNo.H"

kinematicCloud.evolve();
kinematicCloud.info();

fvVectorMatrix UEqn
(
fvm::ddt(rho,U)
+ fvm::div(phi, U)
- fvm::laplacian(mu, U)
==
rho.dimensionedInternalField()*g
+ kinematicCloud.SU(U)
);

solve(UEqn == -fvc::grad(p));

// --- PISO loop

for (int corr=0; corr<nCorr; corr++)
{
volScalarField rUA = 1.0/UEqn.A();

U = rUA*UEqn.H();

phi = fvc::interpolate(rho)*
(
( fvc::interpolate(U) & mesh.Sf())
+ fvc::ddtPhiCorr(rUA, rho, U, phi)
);


adjustPhi(phi, U, p);

for (int nonOrth=0; nonOrth<=nNonOrthCorr; nonOrth++)
{

fvScalarMatrix pEqn
(
fvm::laplacian(rho*rUA, p) == fvc::div(phi)
);

pEqn.setReference(pRefCell, pRefValue);
pEqn.solve();

if (nonOrth == nNonOrthCorr)
{
phi -= pEqn.flux();
}

}

# include "continuityErrs.H"

U -= rUA*fvc::grad(p);
U.correctBoundaryConditions();

}

runTime.write();

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


}

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

return 0;

}

// ************************************************** *********************** //

I think the problem lies with either the source term (quite possibly as this was taken from the coal example) or the way in which the Reynolds number is calculated. I will look into this further.
Thanks,
Paul

The problem was with the value I was specifying for rhoMin in the constant/kinematicCloudProperties file (I had it set to the particle density where it should have been a limiting small value). Fixing this results in the correct Reynolds number being used.

In response to those that have asked for this solver, I am attaching it here along with a simple example case.
I have not had enough time to thoroughly verify that it is working as it should (i.e. it is essentially my first pass at modifying Bernhard Gschaider's icoLagrangianFoam, so that it works in OF2.0), so if anyone uses this and finds any problems, please let me know.

To use the solver extract to your foam applications directory (usually something like /home/username/OpenFOAM/OpenFOAM-2.0.x/applications) and then simply type wmake.

Hi Paul!

I'm intending simulate a multiphase flow in a cilyndrical bioreactor. The tank is composed by a tangential inlet and an internal outlet. I plan to simulate a flow with homogeneous particles with water using one or two-way coupling. Is possible apply icoCoupledKinematicParcelFoam to acess this kind of flow? Is this the newest version of the solver? Do you know any tutorial?

Thank's in advance.

Hi Jonas,

The last time I checked I think icoCoupledKinematicParcelFoam was not imparting the correct amount of momentum on the flow as it should have. So there may still be an issue with the source term and hence with the two way coupling.
I haven't had a chance to reexamine this since (although I hope to look at this again soon).

Your best bet at this stage may be to use Bernhard Gschaider's icoLagrangianFoam (which works well in OF 1.6-ext). it should work for the problem you describe.

Kind Regards,

Paul.