[Amjad Asad]

Hello,
I am trying to couple my solver with Population Balance Equation (PBE) in order to describe the particle size distribution in the flow domain. Unfortunately, the solver is unstable and does not work well.

To caluclate the particle distrubtion, I implement the DQMOM equations ( see equation.pngin OpenFOAM. These equations are taken from this paper done by Silva et al. "Implementation and analysis of numerical solution of the population balance equation in CFD packages"

The implementation in OF used to determine the abscissa, the weight abscissa and the required coffecient can be found in the attachment (pbeDQMOM.H)

{
    for (int acorr=0; acorr< nPBECorr; acorr++)
    {
        forAll(absc[0], cellI)
        {
            // Aggregation (Sa) and breakege (Sb) contribution for cellI
            scalarField Sa = aggK   -> aggregKernel(cellI);
            scalarField Sb = breakK -> breakKernel(cellI); 
           
            // nmom number of moments 
            // nmom=2*np
            // np is the quadarate nodes
                                                
             for(label k=0; k<nmom; k++)
            {
                for(label i=0; i<np; i++)
                {   
                    DQM[k][i]    = (1-k)*pow(absc[i][cellI],k);
                    DQM[k][i+np] = k*pow(absc[i][cellI],k-1);
                }
                source[k] = Sa[k] + Sb[k];
            }

            //- Solve matrix and allocate the results into
            //- source term for DQMOM transport equations
            LUsolve(DQM, source);

            for(label i=0; i<np; i++)
            {
                zeta[i][cellI]  = source[i]; 
                omega[i][cellI] = source[i+np];
            }
        }
              
        //- Solve transport equations for DQMOM
        //- alphac void fraction of contious phase
        surfaceScalarField Alphaf = fvc::interpolate(1.-alphac); 
        
        //beta     = scalar(1.0);
           
        forAll(weight, i)
        {    
            
            surfaceScalarField phir = phid - phi;
            
            //phid disperse phase 
            //phi continous phase 
            //surfaceScalarField phix = fvc::interpolate(1-alphac)*phid    + fvc::interpolate(alphac) *phi;
            
            fvScalarMatrix PBEq1
            (
                fvm::ddt(weight[i])
              + fvm::div(phix,weight[i])
              + fvm::div(phir,weight[i])
              - fvm::div(phir*Alphaf,weight[i]) 
                           ==
                zeta[i]
            );

            fvScalarMatrix PBEq2
            (
                fvm::ddt(eta[i])
              + fvm::div(phix,eta[i])
              + fvm::div(phir,eta[i])
              - fvm::div(phir*Alphaf,eta[i])
              ==             
                omega[i]
            );
            
          
            
            PBEq1.solve();
            PBEq2.solve();
            Info<<"here1" << nl;
            
            absc[i] = eta[i]/(weight[i]);
       
                   
        }
        
        for(label i=0; i<np; i++)
        {
            num += eta[i];
            den += pow(absc[i], 2./3.)*weight[i];
        }
    
        Info<<"here2" << nl;
        ds= CC * num/(den + small);
        Info<<"here3" << nl;
                  
    }        
    
 }

The implementation of aggregation model:

tmp<scalarField> aggMcCoyMadras::aggregKernel (const label cellId) const
{
    scalarField Sa(2*np_, 0.0);
    scalar aux;
    // loop to calculate death and birth due aggregation
    for(label k=0; k < 2*np_; k++)
    {
        for(label i=0; i < np_; i++)
        {
            aux = 0.;
            
            for(label j=0; j < np_; j++)
            {
                aux += ( pow(abscissa_[i][cellId] + abscissa_[j][cellId] , k)
                      - pow(abscissa_[i][cellId] , k) - pow(abscissa_[j][cellId] , k) )
                      * C1.value() * weight_[j][cellId]*weight_[i][cellId];

            }
            
            Sa[k] += aux;
        }
        Sa[k] *= 0.5;
            
     
    }

    return tmp<scalarField>
    ( 
        new scalarField(Sa)
    );
}

The implementation of breakage model:

tmp<scalarField> breakMcCoyMadras::breakKernel (const label cellId) const
{
    scalarField Sb(2*np_, 0.0);
    // loop to calculate death and birth due breakage
    for(label k=0; k < 2*np_; k++)
    {
        for(label i=0; i < np_; i++)
        {
            Sb[k] += 0.5*abscissa_[i][cellId]*phiInf.value()*phiInf.value()
                  * weight_[i][cellId] * pow(abscissa_[i][cellId], k)
                  * ( ni.value()/(k + 1) - 1 );
               
        }
    }

    return tmp<scalarField>
    (
        new scalarField(Sb)
    );
}

the solver:

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 "readGravitationalAcceleration.H"
    #include "createFields.H"
    #include "createPBEVariables.H"
    #include "createFvOptions.H"
    #include "initContinuityErrs.H"

    pimpleControl pimple(mesh);

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

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

    while (runTime.run())
    {
        #include "readTimeControls.H"
        #include "CourantNo.H"
        #include "setDeltaT.H"

        runTime++;
        
 
        
     

        // --- Pressure-velocity PIMPLE corrector loop
        while (pimple.loop())
        {
            #include "UEqn.H"
            #include "pbeDQMOM.H"
              
                
            
            // --- Pressure corrector loop
            while (pimple.correct())
            {
                #include "pEqn.H"
              
            }

            if (pimple.turbCorr())
            {
                turbulence->correct();
                
            }
        }

        runTime.write();

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

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

    return 0;
}

I can compile the implement solver without any errors. But after two or three time steps I get this error massage:

here1
#0  Foam::error::printStack(Foam::Ostream&) at ??:?
#1  Foam::sigFpe::sigHandler(int) at ??:?
#2   in "/lib/x86_64-linux-gnu/libc.so.6"
#3   in "/lib/x86_64-linux-gnu/libm.so.6"
#4  pow in "/lib/x86_64-linux-gnu/libm.so.6"
#5  Foam::pow(Foam::Field<double>&, Foam::UList<double> const&, double const&) at ??:?
#6  void Foam::pow<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>&, Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensioned<double> const&) at ??:?
#7  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::pow<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, Foam::dimensioned<double> const&) at ??:?
#8  Foam::tmp<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> > Foam::pow<Foam::fvPatchField, Foam::volMesh>(Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh> const&, double const&) at ??:?
#9  
 at ??:?
#10  __libc_start_main in "/lib/x86_64-linux-gnu/libc.so.6"
#11  
 at ??:

I think the error comes up because of the height value of the determined absc[i] in the equation den += pow(absc[i], 2./3.)*weight[i] in the code attached here. . Unfortunately I can not recognize the reason of this high value.


It is my pleasure to get your suggestion and advices. Thank you for your help!