Stabilized finite element discretization applied to an operator-splitting method of population balance equations

doi:10.1016/j.apnum.2013.04.001

Applied Numerical Mathematics

Volume 70, August 2013, Pages 58-79

https://doi.org/10.1016/j.apnum.2013.04.001 Get rights and content

Abstract

An operator-splitting method is applied to transform the population balance equation into two subproblems: a transient transport problem with pure advection and a time-dependent convection–diffusion problem. For discretizing the two subproblems the discontinuous Galerkin method and the streamline upwind Petrov–Galerkin method combined with a backward Euler scheme in time are considered. Standard energy arguments lead to error estimates with a lower bound on the time step length. The stabilization vanishes in the time-continuous limit case. For this reason, we follow a new technique proposed by John and Novo for transient convection–diffusion–reaction equations and extend it to the case of population balance equations. We also compare numerically the streamline upwind Petrov–Galerkin method and the local projection stabilization method.

Introduction

Many chemical processes, including polymerization, crystallization, cloud formation, and cell dynamics can be described by population balance equations. Hence, their simulation is required in various applications. A special example is the precipitation process which involves chemical reactions in a flow field. Such processes are modeled by a population balance system [14], consisting of the Navier–Stokes equations for describing the flow field, convection–diffusion problems for describing the chemical reactions, and transport equations for the particle size distribution (PSD). The set of these equations is strongly coupled. Hence, inaccuracies in the concentration of one species directly affects the concentration of all other species. In addition to the coupling of these equations, the main difficulty in simulation is that the PSD depends not only on space and time but also on the properties of particles referred to as internal or property coordinates. Consequently, the dimension of the equation of the PSD is higher than the dimensions of the other equations in the system.

In order to overcome the curse of dimensionality associated with the equation of PSD, we proposed in [1] an operator-splitting scheme. The operator-splitting method reduces the original problem with respect to internal and external directions into a transient transport problem with pure advection and a time-dependent convection–diffusion problem. In applications, most of the problems are convection-dominated and the solution obtained by standard finite element methods exhibit nonphysical oscillations. In this case, stabilization techniques are required in order to get physically sound numerical approximations. In [1], the local projection stabilization (LPS) method has been used to stabilize the space discretization combined with a discontinuous Galerkin (dG) method in the internal coordinate to get the error estimates of the two-step method.

The choice of stabilization parameters plays a critical role in the success of the stabilized methods. The standard techniques to get error estimates for residual based stabilization methods for unsteady convection–diffusion–reaction problems lead to lower bounds on the time step length. The stabilization vanishes in the time-continuous limit. This time step restriction does not appear when stabilization methods like LPS [1], [2] are used. Recently in [12], John and Novo presented a new technique for the case of time-independent coefficients which allows to get the same error estimates without time step restriction. In this paper, we show that this technique can be extended to the case of population balance equations.

Stabilized finite element methods for time-dependent convection–diffusion–reaction problems have been investigated by many authors. The stability property of consistent stabilization methods in the small time step limit have been discussed in [3], [10]. The approach in these studies was to discretize the time-dependent problem in space first and the stabilization is introduced in the semi-discrete problem by using residuals of the time-dependent partial differential equation. Then the stabilization parameters are chosen for the semi-discrete problem. The stabilized semi-discrete problem is then discretized in time by a suitable time-stepping scheme. This procedure results in stabilization parameters that depend only on the mesh width in space since the temporal discretization is performed after the choice of the stabilization parameters. The stability and convergence properties of the SUPG method in space combined with the backward Euler method, the Crank–Nicolson scheme or the second order backward differentiation formula in time for transient transport problem have been studied in [4]. Error bounds in the $L^{2}$ -norm and in the norm of material derivative are obtained under regularity conditions on the data and the stabilization parameters depend only on the mesh size in space. For non-smooth data the bounds are valid if the stabilization parameters are chosen in dependence of the length of the time step. Numerical studies of the different stabilization techniques and a comparison including SUPG method can be found in [7], [17].

On the other hand, if the stabilization parameters are chosen after discretizing the problem in space and time, see [12], [13], then, the stabilization parameters will depend on the time step length. A detailed study of SUPG methods for transient convection–diffusion–reactions equation has been given in [12]. It is shown in the first part of the paper that a finite element discretization in space combined with a backward Euler scheme in time leads to two different choices of stabilization parameters, both depending on the length of the time step. Stability and error estimates are obtained for both choices of stabilization parameters. It is also observed that the stabilization parameters tend to zero as the time step length approaches zero. Furthermore, a special problem where velocity field and reaction do not depend on time has been considered in the second part of the paper. The stabilization parameters are chosen similar to that in the steady-state case. Under certain regularity of the solution, the second part of [12] extended the analysis of [4] and derived error estimates for the $L^{2}$ -norm and the norm of the material derivative with the standard order of convergence. Moreover, an error estimates in the norm of the streamline derivative has been established. The analysis has also been extended to the fully discrete case where backward Euler and Crank–Nicolson methods are used for time discretization. Numerical studies presented in [12], [13] show that this approach leads to solutions which contain nonphysical oscillations for small time steps compared with the approach from [3], [10].

Other results concerning the analysis of stabilized finite element methods for time-dependent convection–diffusion–reaction equations can be found in the literature. We refer to [7], [17] which consider different stabilizing techniques including SUPG. Symmetric stabilization in space combined with the θ-scheme and the second order backward differentiation formula have been considered in [5]. The analysis of discontinuous Galerkin (dG) in time combined with local projection stabilization (LPS) in space has been studied in [2] and in space and time in [9].

The second subproblem in our splitting method is a transport problem with pure advection, so one suitable choice is to approximate it by the discontinuous Galerkin (dG) finite element method [1]. The dG method was first introduced for the neutron transport problem in [19] and then analyzed in [16]. The theoretical analysis of the dG method for scalar hyperbolic equations can be found in [15] and for the space–time dG finite element method in [9]. For an introduction to dG method we refer to [6].

The main focus of the paper consists in deriving error estimates for operator-splitting methods of population balance equations in which the stabilization parameters do not depend on the length of the time step. In particular, we combine the SUPG method in space with the dG method in internal coordinate. For time discretization, a backward Euler time stepping scheme is used. Under certain regularity of the solution, we extend the analysis of [12] to the two-step method of population balance equations and derive error estimates without a lower bound on the time step length. Furthermore, we compare the numerical results with the results obtained by LPS method in space [1].

The remainder of this paper is organized as follows: Section 2 introduces the problem under consideration and the operator-splitting scheme. The SUPG method in space and dG in internal coordinate are introduced in Section 3. In Section 4, the two subproblems are discretized in time using the backward Euler time-stepping scheme. We derive stability estimates for the two-step method and error estimates are given for two different choices of stabilization parameters. Furthermore, error estimates in which the stabilization parameters are independent of the time step length are given. Section 5 presents the combination of LPS in space and dG in internal coordinate. Numerical results illustrating the theory are reported in Section 6 and some conclusions are given in Section 7.

Section snippets

Model problem and operator-splitting scheme

Let $Ω_{x}$ be a domain in $R^{d}$ , $d = 2, 3$ , with polygonal $(d = 2)$ or polyhedral $(d = 3)$ Lipschitz continuous boundary $Γ = \partial Ω_{x}$ , $Ω_{ℓ} = [ℓ_{\min}, ℓ_{\max}] \subset R$ an interval and $T > 0$ the final time. The state of an individual particle in the population balance equation may consists of an external coordinate x, referring to its position in the physical space, and an internal coordinate ℓ, representing the properties of particles such as size, temperature, volume etc. A population balance for a solid process such as

The SUPG and dG method

For a family of admissible and shape regular triangulations $T_{h}$ of the polyhedral domain $Ω_{x}$ , let $V_{h} \subset V$ denote a finite element space of piecewise polynomials of order $r ⩾ 1$ . Note that the inverse inequalities $‖ Δ_{x} v_{h} ‖_{0, K} ⩽ c_{inv} h_{K}^{- 1} ‖ \nabla_{x} v_{h} ‖_{0, K}, ‖ \nabla_{x} v_{h} ‖_{0, K} ⩽ c_{inv} h_{K}^{- 1} ‖ v_{h} ‖_{0, K} \forall v_{h} \in V_{h}$ hold true.

It is well known that the Galerkin method, when applied to problem (5), is unstable and leads to a solution that is globally polluted with huge nonphysical oscillations. A stabilization method becomes necessary. A popular

Time discretization

In this section, we discretize the subproblems in time using the backward Euler time-stepping scheme. Furthermore, we state the main stability bounds for two different choices of stabilization parameters. We consider a uniform partition of the time interval $[0, T]$ with $τ = T / N$ , i.e., $t^{n} = τ n$ , $n = 0, 1, \dots, N$ . After discretizing in time by the backward Euler method, the obtained first order accurate implicit scheme is considered as a two-step method:

First step. For given $u^{n} \in S_{h, k}^{r, q}$ , find ${\tilde{u}}_{h, k}^{n + 1} \in S_{h, k}^{r, q}$

LPS method in space

In this subsection, we explain briefly the LPS method to discretize the subproblems in space. To this end, let $D_{h}$ be the projection space of discontinuous, piecewise polynomials of degree $r - 1$ with $r ⩾ 1$ . Let $D_{h} (K) = {q_{h} |_{K} : q_{h} \in D_{h}}$ be the local projection space and $π_{K} : L^{2} (K) \to D_{h} (K)$ the local $L^{2}$ -projection onto $D_{h} (K)$ . Define the global projection $π_{h} : L^{2} (Ω_{x}) \to D_{h}$ by $(π_{h} v) |_{K} : = π_{K} (v |_{K})$ . The fluctuation operator $κ_{h} : L^{2} (Ω_{x}) \to L^{2} (Ω_{x})$ is given by $κ_{h} : = i d - π_{h}$ where $i d : L^{2} (Ω_{x}) \to L^{2} (Ω_{x})$ is the identity mapping.

We define the

Numerical studies and comparison

This section presents some numerical results for streamline-upwind Petrov–Galerkin, local projection stabilization and discontinuous Galerkin methods applied to population balance equations. All numerical calculations were performed with the finite element package MooNMD [11].

The numerical tests have used $Q_{1}$ , $Q_{2}$ for $V_{h}$ in SUPG method and for $(V_{h}, D_{h})$ the pairs $(Q_{1}^{bubble}, P_{0}^{disc})$ and $(Q_{2}^{bubble}, P_{1}^{disc})$ in LPS method. For discretization in the internal coordinate, discontinuous Galerkin methods of

Conclusion

This paper addressed different ways to obtain error estimates for SUPG and dG finite element methods based on operator-splitting methods applied to population balance equations. For the backward Euler time-stepping scheme, it was shown that standard energy arguments applied to the fully discrete two-step method yield error estimates which couple the choice of stabilization parameters to the length of the time step such that the SUPG stabilization vanishes in the time-continuous limit.

For this

Acknowledgements

The authors acknowledge the financial support from the Deutsche Forschungsgemeinschaft (DFG) under grant MA 4713/2-1, the Federal Ministry of Education and Research (BMBF) under grant 03TOPAA1, Germany, and Kohat University of Science and Technology (KUST-HEC), Pakistan.

References (21)

N. Ahmed et al.
Finite element methods of an operator splitting applied to population balance equations
J. Comput. Appl. Math.
(2011)
N. Ahmed et al.
Discontinuous Galerkin time stepping with local projection stabilization for transient convection-diffusion-reaction problems
Comput. Methods Appl. Mech. Engrg.
(2011)
P.B. Bochev et al.
Stability of the SUPG finite element method for transient advection–diffusion problems
Comput. Methods Appl. Mech. Engrg.
(2004)
E. Burman
Consistent SUPG-method for transient transport problems: stability and convergence
Comput. Methods Appl. Mech. Engrg.
(2010)
E. Burman et al.
Finite element methods with symmetric stabilization for the transient convection–diffusion–reaction equation
Comput. Methods Appl. Mech. Engrg.
(2009)
R. Codina
Comparison of some finite element methods for solving the diffusion–convection–reaction equation
Comput. Methods Appl. Mech. Engrg.
(1998)
M.-C. Hsu et al.
Improving stability of stabilized and multiscale formulations in flow simulations at small time steps
Comput. Methods Appl. Mech. Engrg.
(2010)
V. John et al.
Finite element methods for time-dependent convection–diffusion–reaction equations with small diffusion
Comput. Methods Appl. Mech. Engrg.
(2008)
V. John et al.
Simulations of population balance systems with one internal coordinate using finite element methods
Chem. Eng. Sci.
(2009)
P. Lasaint et al.
On a finite element method for solving the neutron transport equation

There are more references available in the full text version of this article.

Cited by (15)

Crystallization behavior of calcium silicate hydrate in highly alkaline system: Structure and kinetics
2022, Journal of Crystal Growth
During utilization of high silicon solid wastes, a characteristic and complex desilication solution (DSS) generated, remains an important subject to be as precursor of calcium silicate hydrate (C–S–H) crystallization. By using the mixed suspension mixed product removal (MSMPR) crystallizer, the crystallization behavior of C–S–H in DSS was investigated. The effects of temperature, initial Ca/Si ratio, silicate concentrations, solvent alkalinity on the morphology and structure of synthesized C–S–H were investigated through X–ray diffraction, scanning electron microscopy, laser particle size analyzer, and ²⁹Si magic angle spinning nuclear magnetic resonance. The crystallization growth rate and particle size of C–S–H were positively correlated with temperature. The elevation of Ca/Si molar ratio will lead to shortening of the silicate chain to form more monomers. The increase of silicate concentration and solvent alkalinity will result in the formation of shorter silicate chains. This work provides theoretical guidance for precise preparation and design of C–S–H synthesized by a causticization strategy.
Accurate and efficient approximations for generalized population balances incorporating coagulation and fragmentation
2021, Journal of Computational Physics
This study focuses on development of two approaches based on finite volume schemes for solving both one-dimensional and multidimensional nonlinear simultaneous coagulation-fragmentation population balance equations (PBEs). Existing finite volume schemes and sectional methods such as fixed pivot technique and cell average technique have many issues related to accuracy and efficiency. To resolve these challenges, two finite volume schemes are developed and compared with the cell average technique along with the exact solutions. The new schemes have features such as simpler mathematical formulations, easy to code and robust to apply on nonuniform grids. The numerical testing shows that both new finite volume schemes compute the number density functions and their corresponding integral moments with higher precision on a coarse grid by consuming lesser CPU time. In addition, both schemes are extended to approximate generalized simultaneous coagulation-fragmentation problems and retains the numerical accuracy and efficiency. For the higher dimensional PBEs (2D and 3D), the investigation and verification of the numerical schemes is done by deriving new exact integral moments for various combinations of coagulation kernels, selection functions and fragmentation kernels.
Discrete finite volume approach for multidimensional agglomeration population balance equation on unstructured grid
2020, Powder Technology
Citation Excerpt :
Due to complex non-linear nature of the agglomeration population balance equation (PBE), the analytical solutions are only derived for simpler structure kernels such as G∗(t, x, x′) = 1, x + x′ by Gelbard and Seinfeld [16], Fernández-Díaz and Gómez-García [12] and Kaur et al. [23], respectively which leads us to the choice of numerical approximations. Many numerical solutions have been developed by various authors including finite element methods [2,3,5,9,15], Stochastic method [54], method of moments [51–53], least-square spectral method [11], finite volume schemes [13,14,32,36,39,42,44], direct quadrature method of moments [10,25,28,29] or sectional methods like fixed pivot technique (FPT) [49,50] and cell average technique (CAT) [24,35,41,40]. Among all numerical methods, the sectional methods have advantages over the other methods due to its accuracy to predict both various order moments as well as number distribution function [6,34].
In this work, a discrete finite volume scheme (FVS) is developed for approximating a multidimensional agglomeration population balance equation on an unstructured grid. Previous applications of the FVS were restricted only to structured rectangular and triangular grids. The main hindrance of using the unstructured grid is less knowledge of the relation between the quality of solution and the nature of elements. The testing of the proposed FVS is conducted against the FVS on a structured grid using analytically tractable kernels. The numerical results reveal that the FVS with unstructured grid estimates the numerical results with higher accuracy and efficiency than the existing method. The use of unstructured grid offers improved numerical results due to the flexibility of placing the pivot in the space. Moreover, the mixing of components for a bicomponent population is examined using the variance of excess solute for the unstructured grid.
Finite volume approximation of multidimensional aggregation population balance equation on triangular grid
2020, Mathematics and Computers in Simulation
The present work shows the first ever implementation of two-order moments conserving finite volume scheme (FVS) for approximating a multidimensional aggregation population balance equations (PBE’s) on a structured triangular grid. This scheme is based on preservation of the zeroth and conservation of the first order moments. Our main aim is to demonstrate the ability of the FVS to adapt the structured triangular grid well, hence, improves the accuracy of number density function as well as various order moments. The numerical results obtained by the FVS on a triangular grid are compared with the cell average technique. The comparison is also enhanced to illustrate that the FVS with a triangular grid provides the numerical results with higher precision and at lesser computational time as compared to the FVS with a rectangular grid. Additionally, we also study the mixing state of a bicomponent population of clusters (granules) characterized by the normalized variance of excess solute, $χ$ , a parameter that measures the deviation of the composition of each granule from the overall mean. It is shown that the accuracy of the total variance of the excess solute improves when a triangular grid is used in place of a rectangular grid.
Numerical solution of the population balance equation using an efficiently modified cell average technique
2017, Computers and Chemical Engineering
Citation Excerpt :
Another difficulty in solving the PBE numerically is the scarcity of the number of analytically-solved problems that the numerical methods refer to. Most methods are divided into these categories: the pivot-based techniques (Kumar and Ramkrishna, 1996; Kostoglou, 2007; Kumar et al., 2008; Mostafaei et al., 2015), the finite element and finite volume methods (John et al., 2009; Qamar et al., 2009; Ahmeda et al., 2013; Kumar and Kumar, 2013), the Monte Carlo simulation (Lin et al., 2002; Zhao et al., 2007; Xu et al., 2014) and the method of moments (Marchisio and Fox, 2005; Yuan et al., 2012; Bruns and Ezekoye, 2012; Santos et al., 2013). The present method is, in fact, a modified cell average technique which is a fixed-pivot technique (FPT) and a continuation of the previous one (Mostafaei et al., 2015).
In order to increase the accuracy of the modified cell average technique (mCAT) alongside the zeroth and first moments and increase the speed of the technique, a new supplemental solver has been developed in the region of high amount time, a region where the ratio of the zeroth moment to the initial zeroth moment (μ₀(t)/μ₀(0)) is too small or too large depending on the process and the amount of error is high there. The solver is conceptually easy to understand and straightforward for programming. A number of analytically solved problems were simulated by incorporating this solver into the mCAT and the results were compared with the previous solver. Significant improvements for different phenomenon especially the aggregation related systems were observed.
A modified cell average technique for the solution of population balance equation
2015, Journal of Aerosol Science
A new generalized discretization technique, based on the cell average technique (CAT), for solving the population balance equation and predicting the zeroth and first moments has been developed. The proposed algorithm can be easily cast into computer programs, is conceptually straightforward and demonstrated to be computationally economic. Additionally, a new type of mesh, for breakage, aggregation and combined processes is introduced. A new solver, compatible with aforementioned processes, for improving the accuracy of technique has been also proposed. A number of solvable processes were simulated by this technique and compared with the analytical solution and the CAT method. The results were very close to the analytical solution and much more accurate than CAT results.

View all citing articles on Scopus

View full text

Stabilized finite element discretization applied to an operator-splitting method of population balance equations

Abstract

Introduction

Section snippets

Model problem and operator-splitting scheme

The SUPG and dG method

Time discretization

LPS method in space

Numerical studies and comparison

Conclusion

Acknowledgements

J. Comput. Appl. Math.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Comput. Methods Appl. Mech. Engrg.

Chem. Eng. Sci.