Saving specific variables in OpenQBMM in OpenFOAM to enable its use later within OpenFOAM

[kurianjv]

Hey

I was hoping to get some help with respect to saving a variable which is computed in OpenQBMM, like the secondary weights and abscissae at every time step. When I looked at the OpenQBMM 5.0.1 and the issues previously raised on this platform, I found two ways to print variables:

• I can change NO_WRITE to AUTO_WRITE to write the primary/secondary weights/absisscae or sigma in quadratureNode.C. THis writes these variables during runtime at specific intervals that we specify in controlDict. But printing these everytime step and reading it does not sound very efficient.
• Other option was the reconstructPointDistribution.C file which is a postprocessing step to generate NDF at specific points in the domain based on written out moment files.

I cannot figure out how either of these methods can be customized to allow for direct saving variables in OpenFOAM. My thought was to save each of these variables (in a generic way) like under a volScalarField and then use it later in the OpenFOAM solver (rather than within OpenQBMM). This would allow for me to compute NDF at every point in the domain at every time step or calculate other source terms.

Any help would be greatly appreciated!
And thanks in advance

[AlbertoPa]

Hi @kurianjv and welcome!

If I understand your question correctly, you want the entire history of the EQMOM reconstruction (primary and secondary weights, abscissae, sigma). If that is the case, see my comments below.

• I can change NO_WRITE to AUTO_WRITE to write the primary/secondary weights/absisscae or sigma in quadratureNode.C. Tis writes these variables during runtime at specific intervals that we specify in controlDict. But printing these everytime step and reading it does not sound very efficient.

This approach only works for post-processing or for restart purposes. If you need to save all the time steps, it will require to change controlDict so that all the fields that have AUTO_WRITE active to be written, which is going to be inefficient and take a lot of space.

• Other option was the reconstructPointDistribution.C file which is a postprocessing step to generate NDF at specific points in the domain based on written out moment files.

This is not very different from the previous approach. It will require to store all the moments at each time step, and then for each step, invoke the moment inversion to reconstruct the NDF with EQMOM again.

I cannot figure out how either of these methods can be customized to allow for direct saving variables in OpenFOAM. My thought was to save each of these variables (in a generic way) like under a volScalarField and then use it later in the OpenFOAM solver (rather than within OpenQBMM). This would allow for me to compute NDF at every point in the domain at every time step or calculate other source terms.

In OpenFOAM you can force a field to be written at any point in the code with the "write()" method called on the field. For example:

myField.write();

This will create a directory with the time step value as name, and inside it will only store the fields you forced to write. The downside of this approach is that is it not clean: you still get some directories with all the fields that are saved following controlDict. OpenFOAM will not complain (except if you use the incomplete directory as restart point), but Paraview may be less happy.

An alternative for fields registered in the OpenFOAM database can be to use the volFieldValue function object, which you can find documented here with no operation:

• OF.com: https://www.openfoam.com/documentation/guides/latest/doc/guide-fos-field-volFieldValue.html
• OF.org: https://cpp.openfoam.org/v8/classFoam_1_1functionObjects_1_1fieldValues_1_1volFieldValue.html#details

I am not sure about the purpose of this, but if you need to compute source terms based on the EQMOM reconstruction, it may easier on the long run to couple the codes. If you provide some more details on what is the final goal, I can try to guide.

I hope this helps!

[kurianjv]

Hey Alberto

Thank you for the quick reply!
I have previously tried Field.write() and as you mentioned OF printed out the variable but my experience was that OF saves it with a name on how that was computed which made it difficult to visualize on paraview.
I can try to explain myself better, I m trying to couple OpenFOAM and OpenQBMM which is I am trying to access the variables from OpenQBMM in OF.

For example, I was trying to print moment once the populationBalance->solve() is excecuted by
Info << moment.0.populationBalance <<endl;
But there is error that moment is not declared in the solver. But if I declare Moment.0 in createFIelds it reads files that printed out rather than then once computed by PBM.
Another example, I was trying to save the weights from quadratureNode.H to a .H file in OF framework:
scalarField wei = Foam::quadratureNode::weight_;
I get an error saying -
error: a template declaration cannot appear at block scope
template<class scalarType, class vectorType>

I am thinking out loud:)

[AlbertoPa]

Let me try to explain the structure of the code a little bit.

First, quadratureNode is a templated class, so you cannot directly use it as you tried to because you are not specifying the type of field for weights and abscissae. For example, quadratureNode<volScalarField, volScalarField> defines a node at cells centers with scalar weights and abscissa. If you want a node at face centroids, you do quadratureNode<surfaceScalarField, surfaceScalarField>.

This said, if I understand correctly, you just want to access weights and abscissae that are stored already in memory.

The univariatePopulationBalance is selected at runtime and inherits from populationBalanceModel, univariatePDFTransportModel and realizableOdeSolver. The key point here is the inheritance from univariatePDFTransportModel because this object contains the scalarQuadratureApproximation, which is where all the quadrature information is stored.

There is no direct access function to this class from populationBalanceModel because we can leverage the OpenFOAM database to access it. For example, let's say we want to access to the quadrature information in pbeTransportFoam. We can do this as follows:

        const Foam::PDFTransportModels::populationBalanceModels::univariatePopulationBalance& pbe 
            = U.db().lookupObject<Foam::PDFTransportModels::populationBalanceModels::univariatePopulationBalance>
            (
                "populationBalance.populationBalance"
            );

        const auto& quadrature = pbe.quadrature();

        const auto& nodes = quadrature.nodes();

        Info << "Node" << nodes[1].primaryWeight();

To make pbeTransportFoam build with these additions, we need to add

#include "univariatePopulationBalance.H"

to pbeTransportFoam.C

and modify the options file adding the following lines:

    -I$../../../../src/quadratureMethods/populationBalanceModels/lnInclude \
    -I$../../../../src/quadratureMethods/PDFTransportModels/lnInclude \ 
    -I$../../../../src/quadratureMethods/momentSets/lnInclude \ 
    -I$../../../../src/quadratureMethods/realizableOdeSolver \ 
    -I$../../../../src/mappedList \
    -I$../../../../src/mappedPtrList 

The access functions to weights and abscissae are declared here https://github.com/OpenQBMM/OpenQBMM/blob/master/src/quadratureMethods/quadratureNode/quadratureNode/quadratureNodeI.H

Specifically:

• primaryWeight() and primaryAbsicissa() return the primary weight and abscissa in EQMOM and the weight and abscissa in QMOM.
• secondaryWeights() and secondaryAbscisssae() return a PrtList of secondary weights and abscissae for the quadrature node.

Notes

• Access is read-only, you cannot modify the quadrature (and you should not 😄).
• In the code above, I used auto to do this quickly, it would be more readable to use the full type name in a production code.

[kurianjv]

Hey Alberto

Perfect! That worked :)
Thank you so much again.

[AlbertoPa]

You are welcome. I close based on the comment above. Please, reopen at need.

[kurianjv]

Hey Alberto

I was trying to access the source terms used for moment transport equation.
For univariate PDFtrasport model, I understood the OpenQBMM library to

1. Solve implicitly the moment transort equation (only the transient, convection and diffusion terms) inn univariatePDFTransportModel.C
2. Then the source terms, which include nucleation, growth, aggreagation etc, are solved explicitly in the realisableODEsolver for each moment equation

I was trying to save these source terms similarly to volScalarFields in pbetransportFoam (like I was trying to save weights and abscissae which was discussed previously). Could you give some thoughts on how I can access these source terms once they are computed to update the each moment?

[AlbertoPa]

Hi @kurianjv

the source terms are computed on a cell basis, not as fields, because they are solved with an adaptive ODE solver in each cell. The evaluation of the source terms of univariatePopulationBalance is done in

OpenQBMM/src/quadratureMethods/populationBalanceModels/univariatePopulationBalance/univariatePopulationBalance.C

Line 151 in 3f0b4c5

Foam::scalar

However, there is no mechanism in place to access the field of source terms. Why do you need to access it: to change it or to use the value for processing?

[kurianjv]

Hey Alberto

I wanted to get the field for processing, if it is not accessible I will add code to compute it accordingly then. Thank you :)
I had one additional question about the growth model:
With just the Constant growth model activated, OpenQBMM if I am using the QMOM is essentially solving for
$\partial M_kr/\partial t$ + \nabla . (\vec{U}M_k) = \nabla.(D\nabla M_k) + \sum_i G(L_i) k L_i^{k-1} w_i
The source term due to growth is calculated in growthModel.C:
gSourcei = n*Kg(d,lengthbased)sizeorderpow(bAbscissa,sizeOrder-1)

where n is defined as n = node.n(celli, primaryweight, bAbscissa) and d = node.d(celli, bAbscissa). Could you help me in understanding what are n and d ?
And when using constant growth model:
Kg = user defined value (denoted by Cg) when minabscissa <= d <= maxabscissa

Is this understanding correct ?

Thank you so much in advance for taking the time!

[AlbertoPa]

Hi Kurian,

I wanted to get the field for processing, if it is not accessible I will add code to compute it accordingly then. Thank you :)

OK. Just one note: if the source you need is the one for aggregation, recomputing it is very expensive, particularly with EQMOM, compared to the rest of the calculation, because of the number of operations involved. What you could do define extra fields in the OpenQBMM code, with a name, and then access to them using the database, as discussed above. Note that the fields should be properties of some object that persists (like the model itself). You could save one cell value at a time in these fields, and then use the fields later.

I had one additional question about the growth model:
With just the Constant growth model activated, OpenQBMM if I am using the QMOM is essentially solving for
$\partial M_kr/\partial t$ + \nabla . (\vec{U}M_k) = \nabla.(D\nabla M_k) + \sum_i G(L_i) k L_i^{k-1} w_i
The source term due to growth is calculated in growthModel.C:
gSourcei = n*Kg(d,lengthbased)_sizeorder_pow(bAbscissa,sizeOrder-1)

where n is defined as n = node.n(celli, primaryweight, bAbscissa) and d = node.d(celli, bAbscissa). Could you help me in understanding what are n and d ?

n is the number density at that cell:

OpenQBMM/src/quadratureMethods/quadratureNode/quadratureNode/quadratureNode.H

Line 305 in 3f0b4c5

scalar n(const label celli, const scalar& w, const scalar& x) const;

d is the diameter at a cell given the abscissa:

OpenQBMM/src/quadratureMethods/quadratureNode/quadratureNode/quadratureNode.H

Line 295 in 3f0b4c5

scalar d(const label celli, const scalar& x) const;

And when using constant growth model:
Kg = user defined value (denoted by Cg) when minabscissa <= d <= maxabscissa

Is this understanding correct ?

Correct. The pos function is used to define a step function (pos returns a field which is 1 when the argument is positive, zero otherwise).

[kurianjv]

Hey again Alberto

Thank you for the information about files in which both n and d are defined !
My understanding of QMOM is that it writes the source term for growth as sum of quadrature terms:
scalar source = primaryWeight*Kg(d,lengthbased)*sizeOrder*pow(bAbscissa, sizeorder-1)
scalar Sourcefinal = Sourcefinal + source

But in growthModel.C lines 109-139, where the growth term for QMOM is calculated it looks like n is used in place of primaryWeight. When I looked at the quadratureNode/quadratureNode.C to see how n is calculated it seem to be calculated as primaryWeight/pow3(bAbsicssa).

Do let me know I have misunderstood something ?
Also is there reference that you would recommend for understanding the implementation of QMOM in openqbmm? I have been looking at these papers:

• https://www.mdpi.com/2305-7084/2/1/8
• https://www.sciencedirect.com/science/article/pii/S0009250903002112

Thanks again for taking the time!

[AlbertoPa]

Hi Kurian,

you are correct, b is calculated

OpenQBMM/src/quadratureMethods/quadratureNode/quadratureNode/quadratureNode.C

Line 579 in 3f0b4c5

template<class scalarType, class vectorType>

in different ways depending on how the NDF is defined and what the weights represents.

If the distribution is a number density, then the weight is also a number density (number/m^3). However, in some applications this leads to very large numbers, so it is convenient to redefine the equations so that the weight is a volume fraction. This is typically the case when we deal with multiphase flows with large number of bubbles or particles. This is why the n method takes care of returning the actual number density: because the weight may be a number density or a volume fraction or a mass density, depending on the

The growth kernel expects a weight to be a number density, and the abscissa to be a length, that's why the n method is called, so the same kernels can be reused. Note that if volumeFraction = false, then the weight is directly returned:

OpenQBMM/src/quadratureMethods/quadratureNode/quadratureNode/quadratureNode.C

Line 587 in 3f0b4c5

if (!useVolumeFraction_)

About the references, the main ones are listed here: https://www.openqbmm.org/references-and-publications/ In summary:

• The CES paper from Marchisio you listed is the reference we used to implement QMOM in OpenQBMM, but the inversion algorithm we use is from Wheeler because we found it to be more robust.
• The reference from Askari et al. was not used to develop OpenQBMM and came later and the code was modified by the first author to couple it to the multiphase solver. That code did not make it in OpenQBMM because we tend to favor a different approach for bubbly flows (Heylmun et al.), which is what is implemented in OpenQBMM, because it consistently accounts for polycelerity. Coupling QMOM to a two-fluid approach assumes all bubbles in a cell move at the same mean velocity, which is not the case.

For EQMOM (note these two papers were not produced with OpenQBMM: I clarify because someone in a reproducibility study thought they were for some reason. It is OpenQBMM to be based on them):

• T.T. Nguyen, F. Laurent, R.O. Fox, M. Massot, Solution of population balance equations in applications with fine particles: Mathematical modeling and numerical schemes, Journal of Computational Physics. 325 (2016) 129–156. https://doi.org/10.1016/j.jcp.2016.08.017.
• E. Madadi-Kandjani, A. Passalacqua, An extended quadrature-based moment method with log-normal kernel density functions, Chemical Engineering Science. 131 (2015) 323–339. https://doi.org/10.1016/j.ces.2015.04.005.

The two main publications on OpenQBMM are (preprints are on ResearchGate):

• A. Passalacqua, F. Laurent, E. Madadi-Kandjani, J.C. Heylmun, R.O. Fox, An open-source quadrature-based population balance solver for OpenFOAM, Chemical Engineering Science. 176 (2018) 306–318. https://doi.org/10.1016/j.ces.2017.10.043.
• J.C. Heylmun, B. Kong, A. Passalacqua, R.O. Fox, A quadrature-based moment method for polydisperse bubbly flows, Computer Physics Communications. 244 (2019) 187–204. https://doi.org/10.1016/j.cpc.2019.06.005.

A new one is currently under review, on joint size-velocity NDFs. The preprint is here: https://www.researchgate.net/publication/348538579_A_quadrature-based_moment_method_for_the_evolution_of_the_joint_size-velocity_number_density_function_of_a_particle_population