[shuji_type97]

Hi Foamers,

I would like to inquire about the coded source term in OpenFoam. Actually, I'm facing a problem.

When I use "codedSourceTerm" in OpenFoam, there are 3 options "codeAddSup", "codeAddRhoSup", and "codeAddAlphaRhoSup". I am not sure what's the meaning of these options. When I went to check the source code in OpenFoam9(https://cpp.openfoam.org/v9/massSource_8H_source.html), I found the description as below:

// Sources
 
             //- Add a source term to an equation
             FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_SUP);
 
             //- Add a source term to a compressible equation
             FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_SUP);
 
             //- Add a source term to a phase equation
             FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP);

Could anyone please teach me what's the difference between these 3 options and what's the meaning of "equation", "compressible equation", and "phase equation" here?

Meanwhile, I am trying to implement the kinetic model into CFD by source term. I used multiphaseEulerFoam as my solver, liquid and gas phase both are multicomponent. Please find my code for the source term below:

Growth
{
    type            coded;

    selectionMode   all;

    field           biomass.liquid;

    codeInclude
    #{
       #include "fvm.H"
    #};

    codeAddRhoSup
    #{
       Pout<< "**codeAddRhoSup**" << endl;
       
       const vectorField& C = mesh().C();
              
       // Get the required fields
       const volScalarField& B = mesh().lookupObject<volScalarField>("biomass.liquid");// mass fraction of biomass
       const volScalarField& O = mesh().lookupObject<volScalarField>("O2.liquid");// mass fraction of Oxygen
       const volScalarField& S = mesh().lookupObject<volScalarField>("substrate.liquid");// mass fraction of substrate
       const volScalarField& rhoL = mesh().lookupObject<volScalarField>("thermo:rho.liquid");// rho of liquid phase

       // Define the Kinetic parameters
       //dimensionedScalar  mum("mum", dimensionSet(0,0,-1,0,0,0,0), 2.78e-4); // 1/second
       //dimensionedScalar  ks("ks", dimensionSet(1,-3,0,0,0,0,0), 0.05); // kg/m^3
       //dimensionedScalar  ko("ks", dimensionSet(1,-3,0,0,0,0,0), 1e-4); //kg/m^3
       const scalar mum = 2.78e-1; // 2.78e-4
       const scalar ks  = 0;   // 0.05
       const scalar ko  = 0;   // 1e-4

       // Get the timestep
       const scalar deltaT = mesh().time().deltaT().value();
       
       // Creat the growth rate (mu) 
       volScalarField GrowthRate
       (
           IOobject
           (
               "GrowthRate",
               mesh().time().timeName(),
               mesh(),
               IOobject::NO_READ,
               IOobject::AUTO_WRITE
           ),
           mesh(),
           dimensionedScalar("GrowthRate", dimless/dimTime, 0.0)
       );

       // Compute the growth rate
       forAll(C, i)
       {
           GrowthRate[i] = mum*rhoL[i];
       }
       
       // Add to source term
       eqn += GrowthRate*B;
    #};

Here is a simplified kinetic model.

When I use codeAddRhoSup, it works. However, when I use it in codeAddAlphaRhoSup, error is reported:

Could you please teach me why this error happened?

I appreciate your help.

Shuji

[Almanzo2000]

Hi as a fellow coded source user, with how tiny the documentation on coded source is, I'm not really an expert in this case but I can refer you to this post

coded energysource in fvModels

I don't know if this will help you, but as far as I understand codeAddRhoSup is for the compressible solver (e.g. chtMultiRegionFoam) and codeAddSup is for the incompressible one.

best of luck.

[jmilo]

Originally Posted by shuji_type97

Hi Foamers,

I would like to inquire about the coded source term in OpenFoam. Actually, I'm facing a problem.

When I use "codedSourceTerm" in OpenFoam, there are 3 options "codeAddSup", "codeAddRhoSup", and "codeAddAlphaRhoSup". I am not sure what's the meaning of these options. When I went to check the source code in OpenFoam9(https://cpp.openfoam.org/v9/massSource_8H_source.html), I found the description as below:

// Sources
 
             //- Add a source term to an equation
             FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_SUP);
 
             //- Add a source term to a compressible equation
             FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_RHO_SUP);
 
             //- Add a source term to a phase equation
             FOR_ALL_FIELD_TYPES(DEFINE_FV_MODEL_ADD_ALPHA_RHO_SUP);

Could anyone please teach me what's the difference between these 3 options and what's the meaning of "equation", "compressible equation", and "phase equation" here?

Meanwhile, I am trying to implement the kinetic model into CFD by source term. I used multiphaseEulerFoam as my solver, liquid and gas phase both are multicomponent. Please find my code for the source term below:

Growth
{
    type            coded;

    selectionMode   all;

    field           biomass.liquid;

    codeInclude
    #{
       #include "fvm.H"
    #};

    codeAddRhoSup
    #{
       Pout<< "**codeAddRhoSup**" << endl;
       
       const vectorField& C = mesh().C();
              
       // Get the required fields
       const volScalarField& B = mesh().lookupObject<volScalarField>("biomass.liquid");// mass fraction of biomass
       const volScalarField& O = mesh().lookupObject<volScalarField>("O2.liquid");// mass fraction of Oxygen
       const volScalarField& S = mesh().lookupObject<volScalarField>("substrate.liquid");// mass fraction of substrate
       const volScalarField& rhoL = mesh().lookupObject<volScalarField>("thermo:rho.liquid");// rho of liquid phase

       // Define the Kinetic parameters
       //dimensionedScalar  mum("mum", dimensionSet(0,0,-1,0,0,0,0), 2.78e-4); // 1/second
       //dimensionedScalar  ks("ks", dimensionSet(1,-3,0,0,0,0,0), 0.05); // kg/m^3
       //dimensionedScalar  ko("ks", dimensionSet(1,-3,0,0,0,0,0), 1e-4); //kg/m^3
       const scalar mum = 2.78e-1; // 2.78e-4
       const scalar ks  = 0;   // 0.05
       const scalar ko  = 0;   // 1e-4

       // Get the timestep
       const scalar deltaT = mesh().time().deltaT().value();
       
       // Creat the growth rate (mu) 
       volScalarField GrowthRate
       (
           IOobject
           (
               "GrowthRate",
               mesh().time().timeName(),
               mesh(),
               IOobject::NO_READ,
               IOobject::AUTO_WRITE
           ),
           mesh(),
           dimensionedScalar("GrowthRate", dimless/dimTime, 0.0)
       );

       // Compute the growth rate
       forAll(C, i)
       {
           GrowthRate[i] = mum*rhoL[i];
       }
       
       // Add to source term
       eqn += GrowthRate*B;
    #};

Here is a simplified kinetic model.

When I use codeAddRhoSup, it works. However, when I use it in codeAddAlphaRhoSup, error is reported:

Could you please teach me why this error happened?

I appreciate your help.

Shuji

Hey, did you solve the problem? Here is my take on it:
1) in OF9, now you have to be mindful of the kind of problem you're solving.
* for incompressible flows (density doesn't change i.e: flow of liquids in pipes with heat transfer) you probably use "codeAddSup"
* for compressible flows (density changes, i.e: gas reaction) you probably use "codeAddRhoSup" in this case, fvModels adds the density to the source term. You can see this in fvModels.C
* for multiphase flows, you probably use "codeAddAlphaRhoSup" and fvModels will add the density and alpha required on the source term.
2) When you said

I use codeAddRhoSup, it works

I'm not sure it actually worked. Did you see this print out in your output: **codeAddRhoSup**? I bet you didn't and it is because you're not using the right source term for multiphase: "codeAddAlphaRhoSup". Now, when you put your code under "codeAddAlphaRhoSup" you got an error, but I'm fairly certain that is the right way to use it.
3) Now, about this error: the source term has units kg/m3 s. Your "GrowthRate" variable has units of 1/s. This is what the code is complaining about. You might think that because you multiplied by a density it fixes the problems, but apparently "volScalarField" are unitless quantities, so your rhoL doesn't have any units; and OF9 complains about LHS being kg/m3s and the RHS being 1/s. I think the following works:

Growth
{
    type            coded;

    selectionMode   all;

    field           biomass.liquid;

    codeInclude
    #{
       #include "fvm.H"
    #};

    codeAddAlphaRhoSup
    #{
       Pout<< "**codeAddAlphaRhoSup**" << endl;
       
       const vectorField& C = mesh().C();
              
       // Get the required fields
       const volScalarField& B = mesh().lookupObject<volScalarField>("biomass.liquid");// mass fraction of biomass
       const volScalarField& O = mesh().lookupObject<volScalarField>("O2.liquid");// mass fraction of Oxygen
       const volScalarField& S = mesh().lookupObject<volScalarField>("substrate.liquid");// mass fraction of substrate
       const volScalarField& rhoL = mesh().lookupObject<volScalarField>("thermo:rho.liquid");// rho of liquid phase

       // Define the Kinetic parameters
       dimensionedScalar  K_O2("K_O2", dimensionSet(1,-3,0,0,0,0,0), 1.0e-4); //kg/m3
       dimensionedScalar  ub_max("ub_max", dimensionSet(0,0,-1,0,0,0,0), 0.278); //1/s
       dimensionedScalar  rhoB("rhoB", dimensionSet(1,-3,0,0,0,0,0), 1); //kg/m3
       //const scalar mum = 2.78e-1; // 2.78e-4
       //const scalar ks  = 0;   // 0.05
       //const scalar ko  = 0;   // 1e-4

       // Get the timestep
       const scalar deltaT = mesh().time().deltaT().value();
       
       // Creat the growth rate (mu) 
       volScalarField GrowthRate
       (
           IOobject
           (
               "GrowthRate",
               mesh().time().timeName(),
               mesh(),
               IOobject::NO_READ,
               IOobject::AUTO_WRITE
           ),
           mesh(),
           dimensionedScalar("GrowthRate", dimless/dimTime, 0.0)
       );

       // Compute the growth rate
       GrowthRate = ub_max*(O/(O + K_CO2/rho));

       // Add to source term
       eqn += GrowthRate*B*rhoX;
    #};

Specific things to note here:
1) inside the block "codeAddAlphaRhoSup" you can call "rho" and "alpha" and you'll get the density and the volume fraction of the field.
2) for the growth rate I used the Monod model and divided numerator and denominator by the density, leaving it only in terms of the mass fraction of oxygen. GrowthRate has now units of 1/s
3) rhoX is there just to add the missing units. Th RHS of eqn should have units kg/m3s as well.