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

where is the equations for coalChemistryFoam?

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

Like Tree4Likes
  • 1 Post By DustExplosion
  • 2 Post By DustExplosion
  • 1 Post By DustExplosion

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   April 8, 2015, 19:34
Default where is the equations for coalChemistryFoam?
  #1
New Member
 
Join Date: Feb 2015
Posts: 8
Rep Power: 9
musabai is on a distinguished road
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:
Quote:
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.
musabai is offline   Reply With Quote

Old   March 6, 2018, 02:32
Default
  #2
New Member
 
Join Date: Jun 2017
Posts: 12
Rep Power: 7
ms6918 is on a distinguished road
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
ms6918 is offline   Reply With Quote

Old   March 6, 2018, 10:08
Default
  #3
Member
 
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
Posts: 62
Rep Power: 8
DustExplosion is on a distinguished road
You can find some of the source terms in this thesis: http://publications.rwth-aachen.de/r...files/5065.pdf - good luck!
alainislas likes this.
DustExplosion is offline   Reply With Quote

Old   March 7, 2018, 01:59
Default
  #4
New Member
 
Join Date: Jun 2017
Posts: 12
Rep Power: 7
ms6918 is on a distinguished road
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
ms6918 is offline   Reply With Quote

Old   March 8, 2018, 07:05
Default
  #5
Member
 
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
Posts: 62
Rep Power: 8
DustExplosion is on a distinguished road
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
DustExplosion is offline   Reply With Quote

Old   March 8, 2018, 08:44
Default
  #6
New Member
 
Join Date: Jun 2017
Posts: 12
Rep Power: 7
ms6918 is on a distinguished road
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.
ms6918 is offline   Reply With Quote

Old   March 8, 2018, 12:04
Default
  #7
Member
 
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
Posts: 62
Rep Power: 8
DustExplosion is on a distinguished road
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)?
DustExplosion is offline   Reply With Quote

Old   March 8, 2018, 13:38
Default
  #8
New Member
 
Join Date: Jun 2017
Posts: 12
Rep Power: 7
ms6918 is on a distinguished road
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.
ms6918 is offline   Reply With Quote

Old   March 9, 2018, 12:18
Default
  #9
Member
 
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
Posts: 62
Rep Power: 8
DustExplosion is on a distinguished road
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!
DustExplosion is offline   Reply With Quote

Old   March 12, 2018, 03:14
Default
  #10
New Member
 
Join Date: Jun 2017
Posts: 12
Rep Power: 7
ms6918 is on a distinguished road
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
ms6918 is offline   Reply With Quote

Old   March 15, 2018, 13:20
Default
  #11
Member
 
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
Posts: 62
Rep Power: 8
DustExplosion is on a distinguished road
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:

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

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

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

Code:
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.
Attached Images
File Type: png Screenshot from 2018-03-15 15:07:48.png (44.1 KB, 51 views)
File Type: png Screenshot from 2018-03-15 15:13:03.png (3.0 KB, 37 views)
File Type: png Screenshot from 2018-03-15 15:13:38.png (11.7 KB, 35 views)
DustExplosion is offline   Reply With Quote

Old   March 19, 2018, 17:30
Default
  #12
Member
 
Robert
Join Date: Sep 2016
Posts: 32
Rep Power: 8
dussa is on a distinguished road
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:
Code:
//- 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
dussa is offline   Reply With Quote

Old   March 21, 2018, 10:56
Default
  #13
Member
 
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
Posts: 62
Rep Power: 8
DustExplosion is on a distinguished road
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
Attached Images
File Type: png Screenshot from 2018-03-21 12:49:49.png (7.7 KB, 31 views)
File Type: png Screenshot from 2018-03-21 12:50:31.png (12.2 KB, 33 views)
File Type: png Screenshot from 2018-03-21 12:50:44.png (38.7 KB, 32 views)
wbywbywby6 and alainislas like this.
DustExplosion is offline   Reply With Quote

Old   March 21, 2018, 11:01
Default
  #14
Member
 
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
Posts: 62
Rep Power: 8
DustExplosion is on a distinguished road
I just realized that DPMfoam using basicKinematicCollidingCloud does not include thermoCloud/thermoParcels.

The definition in basicKinematicColligingCloud.H is

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

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

Old   March 22, 2018, 13:12
Default
  #15
Member
 
Robert
Join Date: Sep 2016
Posts: 32
Rep Power: 8
dussa is on a distinguished road
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
dussa is offline   Reply With Quote

Old   March 25, 2018, 15:39
Default
  #16
Member
 
Chris Cloney
Join Date: Jun 2016
Location: Halifax, Canada
Posts: 62
Rep Power: 8
DustExplosion is on a distinguished road
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.
alainislas likes this.
DustExplosion is offline   Reply With Quote

Old   March 26, 2018, 18:49
Default
  #17
Member
 
Robert
Join Date: Sep 2016
Posts: 32
Rep Power: 8
dussa is on a distinguished road
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
dussa 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
Guide: Writing Equations in LaTeX on the CFD Online Forums pete Site Help, Feedback & Discussions 27 May 19, 2022 03:19
how to implement k-epsilon equations in ADI method kiamey Main CFD Forum 2 July 17, 2014 17:13
modelling Differential equations in a udf RikardMNorén Fluent UDF and Scheme Programming 2 October 1, 2013 03:36
Riemann invariants of adjoint equations of shallow water equations zqb0929 Main CFD Forum 0 March 15, 2012 00:54
CFD governing equations m.gos Main CFD Forum 0 April 30, 2011 14:21


All times are GMT -4. The time now is 09:01.