Hi all,
this essay is supposed to give a quick overview on numerical methods that are used in chromatography. Please feel free to comment, any suggestions are welcome.
Overview on numerical approaches in chromatography
Introduction
There are various models describing chromatographic processes characterized by different levels of complexities and accounting for different band broadening effects. Considered the most important are the General rate model (GRM), the Lumped rate model with pores (LRMP), the Lumped rate model without pores (LRM) and the ideal model of chromatography, see also Guiochon et al. [8]. These models are described by a system of partial differential equations (PDE) or mixed PDE and algebraic equations.
In order to evaluate the concentration profiles within the column or at the column outlet over time, an explicit closed-form analytical solution to the equations is desirable as it allows for fast and exact computation. However, analytical solutions are only known for specific models and restrictive/simplifying assumptions. The ideal model is the only one of the above mentioned, which is rigorously solved using analytical methods [9]. The method of characteristics can be used to solve the ideal model. Characteristic space-time lines describe constant values over time, depending on the inlet profile. When two characteristic lines intersect, a discontinuity appears which will even happen with continuous injection profiles. The amplitude and progression of the shock layer can be computed based on available mathematical theory [8]. As the ideal model is a huge simplification, its solutions often deviate from what is experimentally observed.
An example for an analytical solution of the more complex GRM was derived by Qamar et al. [12] under the restrictive assumption of linear adsorption isotherms and considering a single component fluid.
Consequently, to provide a comprehensive solver for diverse model assumptions and applications, a numerical approach to approximate the solution is the practice of choice. The basic concept is to discretize the continuous space-time domain of the equations with a grid and solve the equations on the grid points. This leads to approximation errors which typically decay when the grid is refined until the approximation approaches the exact solution. The more spatial grid points (degrees of freedom DOF) and time steps are used, the more computational effort is required. Hence, subjects of interest when comparing different numerical methods are the quality of approximation errors and their behaviour as well as overall computational efficiency.
Numerical schemes – spatial discretization
In numerical simulation the method of lines, i.e., the separation of spatial and temporal discretization has prevailed. A semidiscretization in space is used to derive a system of ordinary differential equations (ODE) which is then discretized in time and algebraically solved by an explicit (e.g. Runge-Kutta) or implicit (e.g. BDF) time integration method.
There are several suitable state-of-the-art methods that can be considered to solve chromatographic models, and the following gives a quick overview on selected methods and some remarks on chromatography-specific considerations.
First, we will discuss different approaches for the spatial semidiscretization.
Finite Differences schemes
Finite Difference methods are based on Taylor’s theorem and use Taylor’s series to derive an arbitrary high order approximation of the spatial derivative. For instance, the spatial first order forward finite difference is derived by
introducing a truncation error of \mathcal{O}(\Delta x). There are numerous ways to reformulate a PDE with finite differences, whereas only some give a numerically stable scheme with rapidly decaying error. A well-known FD scheme is the forward-backward FD by Rouchon et al. [14] which solves the LRM by neglecting the dispersion term in the FD formulation of the LRM and making use of the second order truncation error to approximate the dispersion.
This straightforward approach comes at the cost of a large resulting ODE system due to a relatively fine grid that is needed for accurate approximations. Another drawback is that this method is naturally mass dissipative which can be reduced with grid refinement and a high order FD scheme such as presented by Alhumaize [1]. As mass dissipation is profoundly unphysical, one might consider a mass conserving method especially for large chromatographic systems. There are memetic FD schemes that secure mass conservation but rely on more complex theory and implementation [15]. Another limitation of the FD scheme is the usage of non-rectangular grids in spatial 2-3 D discretization. The computation of the finite differences relies on a stencil of neighbouring grid points (the stencil is defined by the scheme, e.g. the above forward FD only needs the right neighbour while the stencil usually contains more neighbouring grid points). That makes it rather challenging to derive a suitable FD scheme for non-rectangular grids such as curvilinear and unstructured grids where the grid metrics need further consideration. One possible strategy is the transformation into a rectangular grid [16]. A related problem also occurs when the stencil exceeds the domain boundaries, which can be solved by introducing ghost cells outside the domain which is non-trivial especially for complex domains [17].
Finite Volume schemes
For finite volume schemes the grid defines cells ranging between neighbouring grid nodes. Inside these cells the so-called control volumes are defined which give a constant value for each conservative variable inside the cell. This leads to a staircase function defining a local Riemann problem at each cell interface. The solution of the flux function at these interfaces is approximated by a feasible numerical flux of choice, leading to the following semidiscretized formulation (in 1D)
for each control volume u_j , grid spacing \Delta x and the numerical fluxes F at the left and right interface.
This procedure is naturally conservative and monotonicity preserving (see e.g. [2][11]). To enhance the efficiency of the FV scheme, a linear reconstruction of the control volumes is considered for high order FV. This preserves the mass conservation but comes at the loss of the monotonicity property under which every linear high order scheme suffers, known as Godunov’s order barrier theorem [6]. That causes oscillation at discontinuous fronts or extremely steep gradients, known as the Gibbs phenomenon [5]. To overcome this problem, a non-linear mechanism can be built in for the reconstruction such as a slope limiter [2] or a WENO scheme that estimates and weights different polynomials on a neighbouring stencil to estimate a value for the control volume (as derived by von Lieres and Andersson [13] and implemented in CADET). This way, the FV approach allows for the construction of essentially non-oscillatory high order (in this case up to 3) schemes.
Finite Element schemes
Just like with FV, the spatial domain is divided into cells but instead of constant control volumes, finite element schemes define a polynomial of arbitrary high order inside each cell approximating the conservative variables. The great advantage of a polynomial approach within the cells is that FE methods can achieve spectral accuracy with respect to the polynomial degree when the solution is sufficiently smooth [4][10]. Hence a lower number of DOF’s is needed to achieve the same accuracy as a low order (in this case ≤3) scheme, leading to a smaller semidiscretized ODE system to be solved and therefore higher computational efficiency.
There are two major approaches coming with FE schemes. The classical FE is the continuous Galerkin (CG) method where cell interfaces are conditioned to be continuous, leading to a tightly coupled ODE system. The second FE approach is the discontinuous Galerkin (DG) method, which can be thought of as a combination of the FV and the CG scheme. In contrast to CG and in correspondence with FV, the cell interfaces are discontinuous. This way a feasible numerical flux to solve the local Riemann problem can be applied, which adds numerical dispersion to the scheme. This additional artificial dispersion is considered beneficial [3] as a feasible numerical flux yields stabilizing and with refinement decaying dispersion. Compared to that, the CG scheme must be extended by additional dispersive terms (see e.g. a SUPG scheme [20]) to achieve comparable stability properties.
On the other hand, DG results in a larger ODE system than CG since two DOF’s are used at every interface node.
As linear high order schemes, the FE approaches suffer under the Gibbs phenomena at discontinuities. This can either be addressed with for example a shock detection and a low order blending approach. Or it can be neglected when the underlying equations are sufficiently dispersive. Chromatographic models include dispersive terms and multiple band broadening effects that are smoothing the solution. With a dispersive term such as the axial dispersion D_{ax}\frac{\partial^2c}{\partial x^2} in chromatographic models, discontinuities are smoothed instantaneously. Nevertheless, steep gradients might occur e.g. due to a Langmuir isotherm causing self-sharpening fronts [8]. In chromatography we know that another cause of steep gradients can be the injection profile at the column inlet. One can exploit this by applying a highly resolved low order non-oscillatory scheme at the column inlet and a high order DG scheme in the inner column. The DG scheme naturally fits next to a low order FV scheme by exchanging the numerical flux at the interface.
Numerical schemes – temporal integration
The semidiscretization of the underlying equations leads to a system of coupled ODE in time
where R stands for the spatial discretization, the so-called residual.
As mentioned above we can choose an explicit or implicit scheme for the time integration. An explicit scheme relies solely on the current solution and is explicitly evaluated which makes those schemes comparably fast per timestep. The disadvantage of explicit methods is the necessity of a stable time step which can become extremely small for stiff problems (i.e., when the eigenvalues of the linearized semidiscretization are far apart). That problem can dominate the fast evaluation per time step and turn the explicit scheme into a computationally expensive one.
Implicit methods lead to a set of non-linear equations which is solved with Newton’s method requiring the residual Jacobian matrix. In exchange to the larger computational effort per time step, implicit methods allow for large time-step sizes.
In chromatography the semidiscretized ODE system can be very stiff, for example due to a nonlinear Langmuir isotherm during elution causing a stiff front flank [7]. Additionally, the source term at the column inlet has a major influence on the column flow which typically results in a stiff ODE system. Hence an implicit time integration appears advantageous. One might also think of a mixed implicit-explicit scheme to solve stiff and non-stiff parts separately (see e.g. [18][19] for advection-diffusion and diffusion-reaction equations).
It is also of significant importance that the time integration should be of the same order as the spatial discretization. Otherwise, the lower order can dominate the whole method resulting into a lower order scheme with unnecessary high effort for the high order spatial discretization.
References
[1] K. Alhumaizi. Comparison of finite difference methods for the numerical simulation of
reacting flow. Computers & Chemical Engineering, 28(9):1759–1769, 2004.
[2] J. Blazek. Computational Fluid Dynamics: Principles and Applications. Elsevier, Oxford, first edition.
[3] F. Brezzi, B. Cockburn, L. Marini, and E. Süli. Stabilization mechanisms in discontinuous galerkin finite element methods. Computer Methods in Applied Mechanics and Engineering, 195(25):3293–3310, 2006. Discontinuous Galerkin Methods.
[4] B. Cockburn, G. E. Karniadakis, and C.-W. Shu. The Development of Discontinuous
Galerkin Methods. 2000.
[5] J. Foster and F. B. Richards. The gibbs phenomenon for piecewise-linear approximation. The American Mathematical Monthly, 98(1):47–49, 1991.
[6] S. K. Godunov. Finite difference method for numerical computation of discontinuous solutions of the equations of fluid dynamics, (cornell aeronautical lab. transl. from the russian).
Math. Sbornik, Nov. Ser., 47:271–306, 1959.
[7] T. Gu. Modeling and Scale-Up of Size-Exclusion Chromatography, Mathematical Modeling
and Scale-Up of Liquid Chromatography2015, pages 93–104.
[8] G. Guiochon, A. Felinger, and D. G. Shirazi. Fundamentals of Preparative and Nonlinear Chromatography. 2 edition, 2006.
[9] A. S.-M. H. Schmidt-Traub, M. Schulte. Preparative Chromatography.
[10] P. Houston, C. Schwab, and E. Süli. Stabilized hp-finite element methods for first-order hyperbolic problems. SIAM J. Numer. Anal., 37:1618–1643, 2000.
[11] B. Koren. A robust upwind discretization method for advection, diffusion and source terms, pages 117–138. Notes on Numerical Fluid Mechanics. Vieweg, Germany, 1993.
[12] S. Qamar, J. Nawaz Abbasi, S. Javeed, and A. Seidel-Morgenstern. Analytical solutions and moment analysis of general rate model for linear liquid chromatography. Chemical Engineering Science, 107:192–205, 2014.
[13] E. von Lieres and J. Andersson. A fast and accurate solver for the general rate model of column liquid chromatography. Computers & Chemical Engineering, 34(8):1180–1191, 2010.
[14] P. Rouchon, P. Valentin, M. Schonauer, C. Vidal-Madjar, G. Guiochon. Numerical Simulation of band propagation in nonlinear chromatography. J. Phys. Chem. 88 (1985) 2709.
[15] Lipnikov, Konstantin; Manzini, Gianmarco; Shashkov, Mikhail (2014). Mimetic finite difference method. Journal of Computational Physics, 257(), 1163–1227. doi:10.1016/j.jcp.2013.07.031
[16] Mitja Lakner, Igor Plazl. The finite differences method for solving systems on irregular shapes. Computers & Chemical Engineering, Volume 32, Issue 12, 2008, Pages 2891-2896, ISSN 0098-1354, doi: 10.1016/j.compchemeng.2008.02.005
[17] Coco, Armando, Russo, Giovanni (2013). Finite-difference ghost-point multigrid methods on Cartesian grids for elliptic problems in arbitrary domains. Journal of Computational Physics, 241(), 464–501, doi:10.1016/j.jcp.2012.11.047
[18] Uri M. Ascher, Steven J. Ruuth, Raymond J. Spiteri. Implicit-explicit Runge-Kutta methods for time-dependent partial differential equations. Applied Numerical Mathematics, Volume 25, Issues 2–3, 1997, Pages 151-167, ISSN 0168-9274, doi: 10.1016/S0168-9274(97)00056-1.
[19] Verwer, J.G., Sommeijer, B.P. An implicit-explicit Runge–Kutta-Chebyshev scheme for diffusion-reaction equations. SIAM J. Sci. Comput. 2004, 25, 1824–1835.
[20] A. N. Brooks and T. J. R. Hughes, Streamline upwind/Petrov-Galerkin formulations for convection dominated flows with particular emphasis on the incompressible Navier-Stokes equations. Comput. Methods Appl. Mech. Engrg., 32, 199–259 (1982).