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

Why extract Formation Rate term of each species as 'reaction->R(Yi) & Yi'?

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

Like Tree1Likes
  • 1 Post By dlahaye

Reply
 
LinkBack Thread Tools Search this Thread Display Modes
Old   November 23, 2024, 14:20
Default Why extract Formation Rate term of each species as 'reaction->R(Yi) & Yi'?
  #1
New Member
 
Harsh Anand
Join Date: May 2024
Posts: 12
Rep Power: 2
GeekCFD is on a distinguished road
I have successfully extracted the Formation Rate term of each species as 'reaction->R(Yi) & Yi', following the instructions on similar threads. (The full YEqn.H file has been given below.)

I also understand that reaction->R(Yi) returns a fvScalarMatrix object instead of volScalarField object and so cannot be directly used while writing the field values to the volScalarField object.
What I don't understand is that WHY do we need to take the inner product of reaction->R(Yi) with Yi (mass fraction of the i^{th} species) to convert it from fvScalarMatrix object to volScalarField object?

It makes no intuitive sense to take a dot product of 'reaction->R(Yi)' with the mass fraction of the species ('Yi'), as multiplying the formation rate with something between [0,1] will change its value which we don't want.
Someone, please help me understand this complication. I will be grateful for it!


Code:
tmp<fv::convectionScheme<scalar>> mvConvection
(
    fv::convectionScheme<scalar>::New
    (
        mesh,
        fields,
        phi,
        mesh.divScheme("div(phi,Yi_h)")
    )
);

{
    reaction->correct();
    Qdot = reaction->Qdot();
    volScalarField Yt(0.0 * Y[0]);

    forAll(Y, i)
    {
        if (i != inertIndex && composition.active(i))
        {
            volScalarField& Yi = Y[i];

            // Assemble the Yi equation
            fvScalarMatrix YiEqn
            (
                fvm::ddt(rho, Yi)                                  // Unsteady term
              + mvConvection->fvmDiv(phi, Yi)                      // Convection term
              - fvm::laplacian((turbulence->mu() / 1.0 +
                                turbulence->mut() / 0.7), Yi)       // Diffusion term
             ==
                reaction->R(Yi)                                    // Reaction term
              + fvOptions(rho, Yi)                                // Additional source terms
            );

            YiEqn.relax();
            fvOptions.constrain(YiEqn);
            YiEqn.solve(mesh.solver("Yi"));
            fvOptions.correct(Yi);
            Yi.max(0.0);

   if (mesh.time().outputTime())
   {
            // Compute each budget term separately (MODIFICATION START)
            volScalarField unsteadyTerm
            (
               IOobject
               (
                   "unsteadyTerm_" + Yi.name(),
                   mesh.time().timeName(),
                   mesh,
                   IOobject::NO_READ,
                   IOobject::AUTO_WRITE
               ),
               fvc::ddt(rho, Yi)
            );

            volScalarField convectionTerm
            (
               IOobject
               (
                   "convectionTerm_" + Yi.name(),
                   mesh.time().timeName(),
                   mesh,
                   IOobject::NO_READ,
                   IOobject::AUTO_WRITE
               ),
               mvConvection->fvcDiv(phi, Yi)
            );

            volScalarField diffusionTerm
            (
               IOobject
               (
                   "diffusionTerm_" + Yi.name(),
                   mesh.time().timeName(),
                   mesh,
                   IOobject::NO_READ,
                   IOobject::AUTO_WRITE
               ),
               fvc::laplacian((turbulence->mu() / 1.0 +
                               turbulence->mut() / 0.7), Yi)
            );

            volScalarField reactionTerm
            (
               IOobject
               (
                   "reactionTerm_" + Yi.name(),
                   mesh.time().timeName(),
                   mesh,
                   IOobject::NO_READ,
                   IOobject::AUTO_WRITE
               ),
               reaction->R(Yi) & Yi
            );

            // Store terms for visualization
            unsteadyTerm.write();
            convectionTerm.write();
            diffusionTerm.write();
            reactionTerm.write();

            // (MODIFICATION END)
   }

            Yt += Yi;
        }
    }

    Y[inertIndex] = scalar(1) - Yt;
    Y[inertIndex].max(0.0);
}

----------------------------
Best Regards,
GeekCFD
GeekCFD is offline   Reply With Quote

Old   November 28, 2024, 10:20
Default
  #2
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 792
Rep Power: 14
Tobermory will become famous soon enough
It's because fvScalarMatrix is a matrix of coefficients that act on Yi. Think of each diagonal entry in the fvMatrix as the reaction rate divided by the Yi value for that cell.

The & in code:
Code:
            volScalarField reactionTerm
            (
               IOobject
               (
                   "reactionTerm_" + Yi.name(),
                   mesh.time().timeName(),
                   mesh,
                   IOobject::NO_READ,
                   IOobject::AUTO_WRITE
               ),
               reaction->R(Yi) & Yi
            );
is not a dot product. Operator& for a fvScalarMatrix, volScalarField is defined in fvMatrix:
Code:
 template<class Type>
 Foam::tmp<Foam::VolField<Type>>
 Foam::operator&
 (
     const fvMatrix<Type>& M,
     const DimensionedField<Type, volMesh>& psi
 )
 {
     tmp<VolField<Type>> tMphi
     (
         VolField<Type>::New
         (
             "M&" + psi.name(),
             psi.mesh(),
             M.dimensions()/dimVolume,
             extrapolatedCalculatedFvPatchScalarField::typeName
         )
     );
     VolField<Type>& Mphi = tMphi.ref();
  
     // Loop over field components
     if (M.hasDiag())
     {
         for (direction cmpt=0; cmpt<pTraits<Type>::nComponents; cmpt++)
         {
             scalarField psiCmpt(psi.primitiveField().component(cmpt));
             scalarField boundaryDiagCmpt(M.diag());
             M.addBoundaryDiag(boundaryDiagCmpt, cmpt);
             Mphi.primitiveFieldRef().replace(cmpt, -boundaryDiagCmpt*psiCmpt);
         }
     }
     else
     {
         Mphi.primitiveFieldRef() = Zero;
     }
    ...
  
     return tMphi;
 }
i.e. does the cell-wise multiplication of diagonal matrix coeff with the Yi value.
Tobermory is offline   Reply With Quote

Old   November 29, 2024, 05:21
Default
  #3
Senior Member
 
Join Date: Apr 2020
Location: UK
Posts: 792
Rep Power: 14
Tobermory will become famous soon enough
PS if the above is gobbledygook to you (as it was to me, a decade ago), then you'll need to do some research on fvMatrix and how OF solves linear equations. Jasak's thesis (Error analysis & estimation for the FV method) is always a good start. Good luck!
Tobermory is offline   Reply With Quote

Old   November 29, 2024, 07:39
Default
  #4
Senior Member
 
Domenico Lahaye
Join Date: Dec 2013
Posts: 825
Blog Entries: 1
Rep Power: 19
dlahaye is on a distinguished road
Thanks for sharing your insights!

Please have a gobbledygood weekend!
Tobermory likes this.
dlahaye 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
Error in New Surface reaction model (Having multiple reactions) surajkvs OpenFOAM Programming & Development 2 May 23, 2023 22:21
How to get volumetric reaction rate of each species in species transport model Yadu krishnan FLUENT 0 November 28, 2021 00:46
TUI command to access reaction arrhenuis rate Brad_rawl FLUENT 1 September 23, 2021 00:35
Species Reaction Rate RR() kaizhangqmul OpenFOAM 12 September 21, 2021 06:16
UDF changing the rate exponent of a reaction Stefan H Fluent UDF and Scheme Programming 0 September 16, 2009 14:20


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