[musabai]

Hi guys,

I am using coalChemistryFoam. I want to know the equations which coalChemistryFoam uses to calculate the velocities, energy, radiation, etc. Where can I find them? There are a lot of research papers about coal combustion, but the source terms in the governing equations of these papers are quite different. So I want to know what the source term equations in the coal combustion governing equations that coalChemistryFoam uses.

for example, the UEqn.H:

fvVectorMatrix UEqn
(
fvm::ddt(rho, U)
+ fvm::div(phi, U)
+ turbulence->divDevRhoReff(U)
==
rho.dimensionedInternalField()*g
+ coalParcels.SU(U)
+ limestoneParcels.SU(U)
+ fvOptions(rho, U)
);

UEqn.relax();

fvOptions.constrain(UEqn);

if (pimple.momentumPredictor())
{
solve(UEqn == -fvc::grad(p));

fvOptions.correct(U);
K = 0.5*magSqr(U);
}

I don't understand coalParcels.SU(U), what's the equation for this term?

I need your help! Thank you.

[ms6918]

Hello musabai,
Have you got the answer to your question? It would be really helpful if you could tell us where we can get the equations for the source term.
Thank You

[DustExplosion]

You can find some of the source terms in this thesis: http://publications.rwth-aachen.de/r...files/5065.pdf - good luck!

[ms6918]

hello.
Thank You so much for the quick reply. The thesis is very helpful. I have few questions which might sound silly. I am new to Openfoam and working on coalChemistryFoam.
So in the solver, there are some equations which are solved for gaseous phase for ex EEqn.H, pEqn.H etc.
But i could not find where are the equations for particles(cloud). I mean i want to know if there are equations which define the particle state(for ex change in mass of particle). I am really confused.

Thank You again for the reply. Any help is highly appreciated.

Regards

[DustExplosion]

It is pretty complicated but this post here might help: http://www.mydustexplosionresearch.c...ge-flow-chart/

Basically, there is one call out in the solver (i think is cloud_solve or particle_solve or something like that, I do not have my linux computer here). That function is the "start".

Inside there each subsequent class gets checked through in order ReactingMultiphaseCloud/Parcel > ReactingCloud/Parcel > ThermoCloud/Parcel > KinematicCloud/Parcel. If that class/function is found in multiple templates the version from the highest one that it is found in is used (e.g., a function defined in ReactingMultiphaseCloud will be used instead of KinematicCloud).

Lastly, each template particle class has a function called "calc" (e.g., FOAM:ReactingParcel::calc). This is the main solve routine and steps through the "equations". It is a good place to start, but it is important to understand the overall structure, otherwise you might get lost! Each function in calc is defined in one or multiple of the template classes. Hope that helps

[ms6918]

Hello Chris,
Thank You for the quick reply. It really means a lot.
Can you please clarify one doubt. If i want to find the overall source term for any equation, say energy equation, in which file do i have to look as the source term is calculated in both ThermoParcel.H and ThermoCloud.H.

According to my understanding, ThermoParcel.H gives us the source term for a particular parcel and ThermoCloud.H gives us source term for the overall cloud. Am i right?

Kindly correct me.

Again thank you for your valuable reply.

Regards.

[DustExplosion]

From memory, I believe you are correct. Most of the source terms are done in the parcel classes. For example Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer calculates heat transfer to the parcel. I think the cloud classes generally loop through the parcels, but the sources are calculated in the parcel class. One exception might be radiation, where that might be calculated in the cloud class (you would have to check)?

[ms6918]

Thank You for the valuable feedback. Can you please tell me where can i find the equations for lagrangian particles. I can not find the lagrangian case.

Thank you for the constant support.

[DustExplosion]

No problem,

The equations are not computed in one place like the fluid ones. They get calculated across multiple classes and in different functions. You gotta go through them all to figure it out.

For example, the source term from heating is calculated in Foam::scalar Foam::ThermoParcel<ParcelType>::calcHeatTransfer, however the actual equation also depends on the HeatTransferModel (RanzMarshall or NoHeatTransfer).

As another example, the source term for drag is calculated in KinematicParcel (or maybe KinematicCloud?). Part of the equation will be in there, but you will also have to look at the drag model (e.g., SphereDrag, NonSphereDrag, Ext) for the drag coefficient calculation.

Hope that helps!

[ms6918]

Thank You so much for the reply.
I have a small doubt.
In the fils U.Eqn.H, the source term is written as coalParcels.SU(U)
And in KinematicParcel file, 2 source terms are defined as
Su (Explicit momentum source for particle) and dUtrans (Momentum transfer from the particle to the carrier phase)
Can you kindly clarify the difference between the two.
Since it is written that dUtrans is the Momentum transfer from the particle to the carrier phase, shouldn’t this be used in as a source term in U.Eqn.H instead of Su.

Thank you again for the valuable feedback

[DustExplosion]

Alright, I got a chance to dig into the code some more and here are some points that should help.

1) You are correct that dUTrans is the momentum transfer between the particles and fluid. At the bottom of calc this gets added to UTrans. Note that each template class has its own definition of calc. KinematicParcel looks like this:

template<class ParcelType>
template<class TrackData>
void Foam::KinematicParcel<ParcelType>::calc
(
    TrackData& td,
    const scalar dt,
    const label cellI
)
{
    // Define local properties at beginning of time step
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    const scalar np0 = nParticle_;
    const scalar mass0 = mass();

    // Reynolds number
    const scalar Re = this->Re(U_, d_, rhoc_, muc_);


    // Sources
    //~~~~~~~~

    // Explicit momentum source for particle
    vector Su = vector::zero;

    // Linearised momentum source coefficient
    scalar Spu = 0.0;

    // Momentum transfer from the particle to the carrier phase
    vector dUTrans = vector::zero;


    // Motion
    // ~~~~~~

    // Calculate new particle velocity
    this->U_ = calcVelocity(td, dt, cellI, Re, muc_, mass0, Su, dUTrans, Spu);


    // Accumulate carrier phase source terms
    // ~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~~
    if (td.cloud().solution().coupled())
    {
        // Update momentum transfer
        td.cloud().UTrans()[cellI] += np0*dUTrans;

        // Update momentum transfer coefficient
        td.cloud().UCoeff()[cellI] += np0*Spu;
    }
}

I believe this gets added to coalParcels.SU in Foam::KinematicCloud<CloudType>::SU(volVectorField & U) in KinetmaticCloudI.H but you should double check (note the "I" in the file name).

template<class CloudType>
inline Foam::tmp<Foam::fvVectorMatrix>
Foam::KinematicCloud<CloudType>::SU(volVectorField& U) const
{
    if (debug)
    {
        Info<< "UTrans min/max = " << min(UTrans()).value() << ", "
            << max(UTrans()).value() << nl
            << "UCoeff min/max = " << min(UCoeff()).value() << ", "
            << max(UCoeff()).value() << endl;
    }

    if (solution_.coupled())
    {
        if (solution_.semiImplicit("U"))
        {
            const DimensionedField<scalar, volMesh>
                Vdt(mesh_.V()*this->db().time().deltaT());

            return UTrans()/Vdt - fvm::Sp(UCoeff()/Vdt, U) + UCoeff()/Vdt*U;
        }
        else
        {
            tmp<fvVectorMatrix> tfvm(new fvVectorMatrix(U, dimForce));
            fvVectorMatrix& fvm = tfvm();

            fvm.source() = -UTrans()/(this->db().time().deltaT());

            return tfvm;
        }
    }

    return tmp<fvVectorMatrix>(new fvVectorMatrix(U, dimForce));
}

2) You need to be careful with the many definitions of Su in the parcel files. For example Su that is passed into calcVelocity in KinematicParcel.C is zero (it is set as zero three lines before going into the function) as shown above.

3) The real force calculation is computed by Feff = Fcp + Fncp; which includes the coupled and nonCoupled force components. These differ depending on which ParticleForce submodels you are using. In the case of SphereDrag with no gravity you can find the definition of calcCoupled

template<class CloudType>
Foam::forceSuSp Foam::SphereDragForce<CloudType>::calcCoupled
(
    const typename CloudType::parcelType& p,
    const scalar dt,
    const scalar mass,
    const scalar Re,
    const scalar muc
) const
{
    forceSuSp value(vector::zero, 0.0);

    value.Sp() = mass*0.75*muc*CdRe(Re)/(p.rho()*sqr(p.d()));

    return value;
}

Since CalcNonCoupled is not defined, that term is zero when only this submodel is used. And there in calcCoupled you finally have one of your equations (a lot of work to only get one)! This is the inverse of the momentum relaxation timescale.

4) Now it is likely confusing why the momentum relation timescale is calculated and not force on the particle. The reason for this is that the velocity is calculated using an integration scheme in calcVelocity.

const Foam::vector Foam::KinematicParcel<ParcelType>::calcVelocity
(
    TrackData& td,
    const scalar dt,
    const label cellI,
    const scalar Re,
    const scalar mu,
    const scalar mass,
    const vector& Su,
    vector& dUTrans,
    scalar& Spu
) const
{
    typedef typename TrackData::cloudType cloudType;
    typedef typename cloudType::parcelType parcelType;
    typedef typename cloudType::forceType forceType;

    const forceType& forces = td.cloud().forces();

    // Momentum source due to particle forces
    const parcelType& p = static_cast<const parcelType&>(*this);
    const forceSuSp Fcp = forces.calcCoupled(p, dt, mass, Re, mu);
    const forceSuSp Fncp = forces.calcNonCoupled(p, dt, mass, Re, mu);
    const forceSuSp Feff = Fcp + Fncp;
    const scalar massEff = forces.massEff(p, mass);


    // New particle velocity
    //~~~~~~~~~~~~~~~~~~~~~~

    // Update velocity - treat as 3-D
    const vector abp = (Feff.Sp()*Uc_ + (Feff.Su() + Su))/massEff;
    const scalar bp = Feff.Sp()/massEff;

    Spu = dt*Feff.Sp();

    IntegrationScheme<vector>::integrationResult Ures =
        td.cloud().UIntegrator().integrate(U_, dt, abp, bp);

    vector Unew = Ures.value();

    // note: Feff.Sp() and Fc.Sp() must be the same
    dUTrans += dt*(Feff.Sp()*(Ures.average() - Uc_) - Fcp.Su());

    // Apply correction to velocity and dUTrans for reduced-D cases
    const polyMesh& mesh = td.cloud().pMesh();
    meshTools::constrainDirection(mesh, mesh.solutionD(), Unew);
    meshTools::constrainDirection(mesh, mesh.solutionD(), dUTrans);

    return Unew;
}

Following through this code and remembering that massEff = mass of the parcel (unless effective mass is being accounted for) and Feff.Su = 0 should give you an equation that look something like 3.60 to 3.62 in the attached. Note that 3.60 can be derived from 3.50, 3.57, and 3.58 which are also attached as images.

[dussa]

Hi Chris and ms6918 (and everyone else who is reading),

I am also working with lagrangian solvers (although I am looking at DPMFoam and not the coalChemistryFoam solver). Reading through your discussions at appears though that you may have some insights that would help in my work. I am looking to add internal heating to particles in the DPMFoam solver. I have seen how heat transfer can be implemented with the use of the RanzMarshall model, however in my particular case I need to define an internal heating value for the particles. The internal energy from the particles will then be transferred to the surrounding liquid phase which is moving past the pebbles.

My question is; where can I actually define an energy source term for pebbles? I have looked through ThermoParcel.H (after seeing your example for the location of the heatTransfer function). It seems that in this file there is a definition for explicit particle enthalpy:

//- Calculate new particle temperature
        template<class TrackData>
        scalar calcHeatTransfer
        (
            TrackData& td,
            const scalar dt,           // timestep
            const label celli,         // owner cell
            const scalar Re,           // Reynolds number
            const scalar Pr,           // Prandtl number - surface
            const scalar kappa,        // Thermal conductivity - surface
            const scalar NCpW,         // Sum of N*Cp*W of emission species
            const scalar Sh,           // explicit particle enthalpy source
            scalar& dhsTrans,          // sensible enthalpy transfer to carrier
            scalar& Sph                // linearised heat transfer coefficient
        );

I wonder if you have had any experience with this term, and if you think this is where I should start if I want to define an internal energy for the particles in my system? Perhaps there is already a simple way of including particle heating for particle types that I am unaware of in the existing solver framework?

I am assuming I will need to create an altered DPMFoam solver which includes additional thermoParcel/Cloud src code, but at the moment I am lost as to how to proceed. If you have any suggestions I would greatly appreciate it.

Kind regards,
Robert

[DustExplosion]

Hi Dussa,

You are in the correct spot, but you need to make sure you have the correct variables. The variable Sh is an input not an output. If you are using ThermoParcels this value is set to 0.0 right before calling calcHeatTransfer.

The equations for particle heating are attached. Note that Qv,p and Qs,p are zero unless devolatilization and surface reaction are included.

For compiling, it should follow pretty closely to this group that changed the force routines: http://www.tfd.chalmers.se/~hani/kur...aanMarkale.pdf

[DustExplosion]

I just realized that DPMfoam using basicKinematicCollidingCloud does not include thermoCloud/thermoParcels.

The definition in basicKinematicColligingCloud.H is

#ifndef basicKinematicCollidingCloud_H
#define basicKinematicCollidingCloud_H

#include "Cloud.H"
#include "KinematicCloud.H"
#include "CollidingCloud.H"
#include "basicKinematicCollidingParcel.H"

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

namespace Foam
{
    typedef CollidingCloud
    <
        KinematicCloud
        <
            Cloud
            <
                basicKinematicCollidingParcel
            >
        >
    > basicKinematicCollidingCloud;
}

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

#endif

This is compared to the following definition for coalCloud

#ifndef coalCloud_H
#define coalCloud_H

#include "Cloud.H"
#include "KinematicCloud.H"
#include "ThermoCloud.H"
#include "ReactingCloud.H"
#include "ReactingMultiphaseCloud.H"
#include "coalParcel.H"

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

namespace Foam
{
    typedef ReactingMultiphaseCloud
    <
        ReactingCloud
        <
            ThermoCloud
            <
                KinematicCloud
                <
                    Cloud
                    <
                        coalParcel
                    >
                >
            >
        >
    > coalCloud;
}

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

#endif

Notice that ThermoCloud is missing.

I think your particles are non-heating (it is assumed that they are in equilibrium with the fluid). You will need to add thermoCloud to the definition to create a new particle class or use a different class that has it included.

[dussa]

Hi Chris,

Thankyou for your quick reply!
I noticed that the Sh term was set to zero, but I was not really sure why this was the case. I was hoping that it was going to be the term that I could define as the particle internal heating value, but I suppose I may have to come up with a different method.

Yes that is the case at the moment with DPMFoam, the particles are at the same temperature as the fluid. However I have seen that there is still the option to include a heat transfer model (like RanzMarshall) in the kinematicCloudProperties file, so perhaps there is still partial heat treatment in the DPMFoam model.

Either way, I think after reading your last message that the way to include heating in the DPMFoam solver is to create my own version of the solver and include the ThermoCloud libraries in a similar way to them being included in coalCloud. Seeing as there is still no term for internal heat generation in the coalChemistryFoam (from what I can see) I might have to alter the ThermoCloud libraries to include some energy source term for the particles if that functionality is missing all together. Do you see any issues with this?

I saw your website with the flow diagram of coalChemistryFoam, it is a great way to visually see the dependencies of the different portions of the code and is very helpful.

Kind regards,
Robert

[DustExplosion]

Thanks Robert,

I think the Sh is there for extension to other solvers that may modify the energy through other mechanisms (e.g., reaction).

Your process sounds correct for creating the new solver. One approach might be to rewrite equation 3.63 and 3.64 in my previous post to their form with a non-zero Biot number. Then you would modify calcHeatTransfer to follow your new derivation.

You could also do it as a separate source term as well. You just need to make sure the heat transfer to the fluid is adjusted accordingly.

[dussa]

Thanks for your advice, I will keep on going along with this methodology then.

Also, I think you are correct with the Sh term, it does seem to be predominantly called in the reacting models such as ReactingParcel and ReactingMultiphaseParcel, although as it is already in the ThermoParcel as well, perhaps if I alter its use in my own version of the lagrangian libraries it could be used as the definition of the internal heating power of my particles.

That was my initial thought, to use a separate source term that is read from a dictionary (and eventually from a file for each time folder, as eventually I will be aiming to read the power of each individual pebble which will be generated with another code).

Thankyou for your help, I will probably start more threads as I get deeper into the issue.
Kind Regards,
Robert