Handle the negative growth rate and limited abscissa

[xinzhenQin]

I based on OF 2106 and chose the corresponding version of QBMM.
First of all, when calculating particle growth, I found that when the growth rate is negative, there is a problem in calculating ODE.
Second, for a positive growth rate, when maxAbscissa is set, the ODE solver will be stuck near maxAbscissa. Even the pureGrowth calculation example that comes with pbeFoam has this problem

[AlbertoPa]

Thank you for reporting, @xinzhenQin

In general, negative growth is not supported out of the box. It requires removing the quadrature abscissae that approach zero and using the remaining nodes only, which was not implemented in the current algorithm. I marked this part as a feature request.

I will investigate the second part, concerning positive growth and let you know.

• Investigate/correct issue with limited abscissa leading to infinite loop in realizable ODE solver
• Split out growth integration
• Manage negative growth

[AlbertoPa]

I confirm the problem: when setting a maximum value for the abscissa, the source term becomes zero, as expected once such value is reached. In the case of pure growth, this leads to no change in the moments, which causes the ODE solver to loop indefinitely.

I tested a correctio against the pureGrowth case with maxAbscissa set to a sufficiently small value to trigger the situation, and it now runs. Cases for pbeFoam are also correctly executed and validated.

See commit 63c6677, which works with OpenFOAM v2212.

Managing negative growth will take some more time.

[xinzhenQin]

Thank you very much for your quick response！But I still wonder if I can make the corresponding changes in the QBMM version corresponding to of2016.

[AlbertoPa]

Yes, you could make that change: see line 321 and following here: 63c6677.
However, we only maintain the version for the current OpenFOAM release, so these changes won't be backported.

[xinzhenQin]

Thanks, I have run successfully on of2106 against the modified version. Thank you very much for your work. But I also found that sometimes the ODE solver reports problems that are not guaranteed to be stable.

[AlbertoPa]

Could you please provide more details on the issues you have found? If you could point to a small case to reproduce them, it would be great (also, not that for negative growth, I have not pushed a solution yet, and how the problem is addressed depends on if QMOM or EQMOM is used. In the first case, you would remove quadrature nodes when the abscissa is null, in the second case the strategy in the paper from Nguyen et al. (2016) can be used).

[xinzhenQin]

Ok, I'll put together a description of the problem and send it up here later today @AlbertoPa

[xinzhenQin]

I'm sorry, I think I'm having some architectural issues. For example, I couple the univariate populationBlanceModel with a solver, but the information of each order moiments such as volScalarField--"momnet.o.populationBalance" are stored under the univariatePDFTransportModel class. Since the populationBlanceModel class does not contain this information, I cannot obtain this field information through populationBlanceModel object created in my solver. I think if this problem can be solved, QBMM can be better coupled with CFD！

[AlbertoPa]

@xinzhenQin this is not the case. Fields of moments and nodes are accessible via the OpenFOAM database. Once the simulation is running, you can retrieve them as normally done with any field. Both are contained in the quadratureApproximation object, so retrieve the population balance model object as shown here

#166 (comment)

then retrieve the quadratureApproximation object, and from there you have access to all the data you need.

Access functions available from the quadratureApproximation object can be seen here: https://github.com/OpenQBMM/OpenQBMM/blob/master/src/quadratureMethods/quadratureApproximations/quadratureApproximationI.H

In general, OpenQBMM heavily relies on OpenFOAM structures and facilities to avoid duplicating code and simplify maintenance, so what can be done in an OpenFOAM application in terms of data access is typically doable with OpenQBMM classes. I would suggest you take some time to explore the code if you plan to use OpenQBMM and couple it to custom solvers.

[xinzhenQin]

thank you very much @AlbertoPa
But this approach relies on creating a concrete PBE model object in the solver, which seems to contradict the mechanism chosen at runtime, although it is feasible for specific problems. This is very similar to the turbulence model, but the abstract classes of the turbulence model provide an interface to rich information. Perhaps it would be more convenient for the abstract class of PBE to also provide corresponding functions.

[AlbertoPa]

Currently, populationBalanceModel is just a wrapper class (not really aware of what quadrature type the actual pbe is going to use) to allow the runtime selection of different types of population balance models. However, the actual calculations happen in the univariatePDFTransportModel and in the specialized version of it univariatePopulationBalanceModel.

The logic is that the PDFTransportModel is associated to a quadratureApproximation, consequently it returns the corresponding quadrature and moments.

[xinzhenQin]

Thank you for your work！@AlbertoPa

[xinzhenQin]

Since we can't currently handle the problem of abscissa less than 0 caused by negative growth rates, I chose to add a minAvscissa to avoid abscissa less than 0. When I set up a pure growth study as follows, a new error was raised：

1. Modify the initial condition: moment.0.populationBalance--"3.125e12"; moment.1.populationBalance--"6.25e6"; moment.2.populationBalance--"12.5"; moment.3.populationBalance--"2.5e-5"; moment.4.populationBalance--"5e-11";
2. Modify the populationaBalanceProperties file:
   growth on;
   growthModel
   {
   growthModel constant;
   minAbscissa 1.0e-9;
   Cg Cg [0 0 -1 0 0 0 0 ] -9;
   }
   error:
   Solving source terms in realizable ODE solver.
   Not realizable, adjusting local timestep.
   This may take a while.

--> FOAM FATAL ERROR: (openfoam-2106)
Reached minimum local step in realizable ODE
solver. Cannot ensure realizability.

From void Foam::realizableOdeSolver<momentType, nodeType>::solve(Foam::realizableOdeSolver<momentType, nodeType>::quadratureType&, Foam::label) [with momentType = Foam::moment<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::quadratureNode<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> > >; nodeType = Foam::quadratureNode<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> >; Foam::realizableOdeSolver<momentType, nodeType>::quadratureType = Foam::quadratureApproximation<Foam::moment<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::quadratureNode<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> > >, Foam::quadratureNode<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::GeometricField<Foam::Vector<double>, Foam::fvPatchField, Foam::volMesh> > >; Foam::label = long int]
in file /home/user/QXZ/applications/OpenQBMM/OpenQBMM-OpenQBMM_7.0.0/src/quadratureMethods/realizableOdeSolver/realizableOdeSolver.C at line 276.

FOAM aborting

#0 Foam::error::printStack(Foam::Ostream&) at ??:?
#1 Foam::error::exitOrAbort(int, bool) at ??:?
#2 Foam::realizableOdeSolver<Foam::moment<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::quadratureNode<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::GeometricField<Foam::Vector, Foam::fvPatchField, Foam::volMesh> > >, Foam::quadratureNode<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::GeometricField<Foam::Vector, Foam::fvPatchField, Foam::volMesh> > >::solve(Foam::quadratureApproximation<Foam::moment<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::quadratureNode<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::GeometricField<Foam::Vector, Foam::fvPatchField, Foam::volMesh> > >, Foam::quadratureNode<Foam::GeometricField<double, Foam::fvPatchField, Foam::volMesh>, Foam::GeometricField<Foam::Vector, Foam::fvPatchField, Foam::volMesh> > >&, long) at ??:?
#3 Foam::PDFTransportModels::univariatePDFTransportModel::solve() at ??:?
#4 ? at ??:?
#5 __libc_start_main in /lib64/libc.so.6
#6 ? at ??:?
Aborted (core dumped)

[AlbertoPa]

You cannot simply bind the abscissa so it does not become negative if you have a negative growth term. The zero-order moment at some point will become inconsistent with the other moments, which changed due to negative growth, while m0 did not. In QMOM, you must take care of removing the quadrature nodes whose abscissa is exactly zero. In EQMOM you can implement the approach discussed in:

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.

which refers to:

M. Massot, F. Laurent, D. Kah, S. De Chaisemartin, A Robust Moment Method for Evaluation of the Disappearance Rate of Evaporating Sprays*, SIAM Journal on Applied Mathematics. 70 (2010) 3203–3234.