Efficient Simulation of Chromatographic Processes Using the Conservation Element/Solution Element Method

Abstract

Chromatographic separation processes need efficient simulation methods, especially for nonlinear adsorption isotherms such as the Langmuir isotherms which imply the formation of concentration shocks. The focus of this paper is on the space–time conservation element/solution element (CE/SE) method. This is an explicit method for the solution of systems of partial differential equations. Numerical stability of this method is guaranteed when the Courant–Friedrichs–Lewy condition is satisfied. To investigate the accuracy and efficiency of this method, it is compared with the classical cell model, which corresponds to a first-order finite volume discretization using a method of lines approach (MOL). The evaluation is done for different models, including the ideal equilibrium model and a mass transfer model for different adsorption isotherms—including linear and nonlinear Langmuir isotherms—and for different chromatographic processes from single-column operation to more sophisticated simulated moving bed (SMB) processes for the separation of binary and ternary mixtures. The results clearly show that CE/SE outperforms MOL in terms of computational times for all considered cases, ranging from 11-fold for the case with linear isotherm to 350-fold for the most complicated case with ternary center-cut eight-zone SMB with Langmuir isotherms, and it could be successfully applied for the optimization and control studies of such processes.

1. Introduction

Chromatographic separation processes are used for the separation of temperature-sensitive mixtures and mixtures of components with very similar physical properties, making them difficult to separate via other cheaper methods. These types of processes are frequently used in the chemical and the pharmaceutical industries ranging from small-scale batch separations of highly valuable pharmaceutically active compounds to large-scale continuous separations of isomer mixtures in the petroleum industries [1].

Mathematical modeling of such processes usually leads to a system of partial differential equations (PDEs) [1]. These mathematical models are used for the design, optimization, and control of chromatographic processes. Depending on the modeling assumptions, different types of models are available and frequently used including equilibrium models, which assume thermodynamic equilibrium between the solid and the fluid phase, or mass transfer models with finite mass transfer resistance between both phases.

For negligible axial dispersion, the equilibrium model admits an analytical solution for piecewise constant initial and boundary conditions and for certain types of linear and nonlinear adsorption isotherms including the well-known Langmuir isotherm. This analytical approach is the core of the so-called equilibrium theory [2,3], which has become an important tool for the conceptual design of chromatographic processes [4,5,6,7,8]. For the mass transfer models, analytical solutions are only possible under further restrictions [9,10,11]. Thus, in both cases, the numerical solution of the underlying PDE system is an important tool for the validation of the analytical results and to overcome the limitations posed by the analytical treatment.

In general, the numerical simulation of chromatographic columns is challenging due to the occurrence of steep concentration fronts, which was the subject of many papers [12,13,14]. For the real-time control and the optimization of the chromatographic processes, fast and reliable methods are required. In practice, the method of lines (MOL) [15] is often used for the solution of the PDE system using finite differences (FDs), finite elements (FEM), or finite volumes (FVM) for the spatial discretization. The simplest choice is a first-order finite volume discretization leading to the popular cell models (CELL), which can then be solved with some standard integrators. The efficiency and accuracy can be further enhanced, with more advanced high-resolution schemes and adaptive grids [13,14,16]. For a recent review of the available numerical methods and simulation tools with many references, the reader is referred to Chapter 6.4 of [1].

An interesting alternative is the so-called space–time conservation element/solution element (CE/SE) method which was originally developed by Chang [17,18,19] for similar problems in fluid dynamics. Later, the CE/SE method was also successfully applied to chromatographic problems by Jørgensen and coworkers [12,20]. Focus was on a mass transfer model for binary simulated moving bed (SMB) processes with linear and nonlinear adsorption isotherms. It was shown [12] that the solutions obtained using the CE/SE method are comparable in accuracy with a MOL approach using a fifth-order high-resolution weighted essentially non-oscillatory (WENO) discretization scheme but need a significantly shorter calculation time. More recently, results were extended in [21] by combining the CE/SE method with a continuous prediction technique. Again, focus was on the mass transfer model for binary SMB processes with linear and nonlinear adsorption isotherms, and again it was shown for a benchmark problem that the CE/SE method is very fast compared to an MOL and a wavelet collocation approach. This motivated us to investigate possible applications of the CE/SE method to the other important class of mathematical models of chromatographic columns, i.e., equilibrium models, and to extend its field of application to a challenging class of ternary SMB processes, i.e., so-called center-cut separations, where the process objective is to isolate a component with intermediate affinity to the solid phase from a mixture of at least three components. This type of separation process plays an important role in the separation of natural products and in biotechnology [22]. In addition, to the best of our knowledge, the first dynamic simulation of a substantially nonlinear eight-zone SMB process is presented, which is much more challenging than the linear [23] or the weakly nonlinear [24] cases which were previously considered [23].

The outline of the paper is as follows: first, the different mathematical models of chromatographic process are presented in Section 2, and the CE/SE method is described in Section 3. Subsequently, numerical studies are shown in which the cell model (MOL) and CE/SE method are compared in terms of efficiency for different process configurations, models, and isotherms in Section 4; lastly, discussion and conclusions are presented in Section 5.

2. Mathematical Models of Chromatographic Processes

The starting point is the linear driving force (LDF) model of a single chromatographic column, which assumes a finite mass transfer resistance between the solid and the fluid phase. In the present paper, it is used for the modeling of single-column batch processes and multi-column continuous simulated moving bed (SMB) processes. It consists of three equations—one partial differential equation (PDE), one ordinary differential equation (ODE), and one algebraic equation (AE) for each of the components and each of the columns.

$$\{\begin{array}{c}\begin{array}{c}\epsilon \frac{\partial C_{i,k}}{\partial t}+\left(1-\epsilon \right)\frac{\partial q_{i,k}}{\partial t}+\epsilon v_k\frac{\partial C_{i,k}}{\partial z}=\epsilon D_{ax}\frac{{\partial }^2{C}_{i,k}}{\partial z^2}\\ \frac{\partial q_{i,k}}{\partial t}=k_m{}_{i,k}\left(q_{i,k}^{*}\left(\mathit{C}\right)-q_{i,k}\right)\end{array} \\ q_{i,k}^{*}=f\left(\mathit{C}\right)\end{array},$$

(1)

where ε is the column void fraction, v is the liquid phase velocity (m·s⁻¹), D_ax is the axial dispersion coefficient (m²·s⁻¹), k_m is the mass transfer coefficient between the two phases (s⁻¹), C is the liquid-phase concentration (g·L⁻¹), q is the solid-phase concentration (g·L⁻¹), q* is the solid-phase concentration at the interphase boundary in equilibrium with the liquid phase (g·L⁻¹), t and z are time and spatial coordinates, respectively, (s) and (m), and i and k are the component and column indices, respectively. This model formulation assumes isothermal operation, a constant void fraction, and a constant mobile-phase velocity inside each of the columns. The apparent axial dispersion coefficient D_ax lumps together all effects leading to band broadening in addition to the finite mass transfer resistance, which is known to have a similar effect [2]. Below, only the limiting case when D_ax $\to$ 0 is considered which is valid for efficient columns with a high number of theoretical stages N. The algebraic equation describes the thermodynamic equilibrium between the solid and the liquid phase and represents the adsorption isotherm. In the present work, two types of adsorption isotherms are used:

$$q_{i, k}^{*}=H_{i,k}{C}_{i,k},$$

(2)

$$q_{i, k}^{*}=\frac{H_{i,k}{C}_{i,k}}{1+{{\displaystyle \sum }}_{r=1}^{N_{\mathrm{comp}}}{b}_{r,k}{C}_{r,k}},$$

(3)

where the former is linear, and the latter is the well-known Langmuir isotherm. H denotes the adsorption Henry coefficient, b is the retention factor (L·g⁻¹), and r is the component index. In the limiting case where the mass transfer is instantaneous (i.e., k_m,i $\to \infty$ and $q_{i,k}\to q_{i,k}^{*}$) and there is negligible dispersion, the ideal equilibrium model [1] is obtained.

$$\epsilon \frac{\partial C_{i,k}}{\partial t}+\left(1-\epsilon \right)\frac{\partial q_{i,k}^{*}\left(\mathit{C}\right)}{\partial t}+\epsilon v_k\frac{\partial C_{i,k}}{\partial z}=0.$$

(4)

The loading of an empty column is now considered as a test scenario. For the present isotherms, this scenario is more challenging than the regeneration step due to the occurrence of steep concentration fronts and is, therefore, used as a benchmark problem. Frequently applied pulse injections can be considered as a combination of loading and regeneration steps and are, therefore, also included to some extent. The initial conditions (ICs) are zero, i.e., we start every simulation with empty columns.

$$\{\begin{array}{c}{C}_{i,k}\left(z, t=0\right)=0\\ q_{i,k}\left(z, t=0\right)=0\end{array}.$$

(5)

Boundary conditions (BCs) are of the Dirichlet type,

$$C_{i,k}\left(z=0 , t\right)=C_{z=0^{-}},$$

(6)

where $C_{z=0^{-}}$ is the liquid-phase concentration just before the column inlet. For the SMB processes, these BCs follow from the material balances of the in- and outlet ports as later discussed in detail.

3. Conservation Element/Solution Element (CE/SE) Method

The most distinctive feature of the CE/SE method developed by Chang [17,18,19] is that time and space are treated in the same manner through the so-called conservation element and solution element. It is based on the integral formulation of the conservation laws and can resolve the typical steep concentration fronts of nonlinear chromatography. It is an explicit time-marching method for the solution of PDEs, i.e., the value of the quantity of interest or so-called state variable (liquid-phase C and solid-phase q concentrations in chromatography) in the current time instant and its spatial derivative only depend on their values at the previous time instants (Figure 1a).

The model in Equation (1) with D_ax $\to$ 0 can be presented in the following vector form (the indices are omitted for simplicity):

$$\mathit{u}=\left(\begin{array}{c}C\\ q\end{array}\right),$$

(7)

$$\mathit{f}=\left(\begin{array}{c}vC\\ 0\end{array}\right),$$

(8)

$$\mathit{p}=\left(\begin{array}{c}-\frac{1-\epsilon }{\epsilon }{k}_m\left(q^{*}\left(\mathit{C}\right)-q\right)\\ k_m\left(q^{*}\left(\mathit{C}\right)-q\right)\end{array}\right),$$

(9)

where u is vector of state variables, f is vector of fluxes, and p is vector of source terms. Equation (1) can then be written in the following form:

$$\frac{\mathit{\partial }\mathit{u}}{\mathit{\partial }\mathit{t}}+\frac{\mathit{\partial }\left(\mathit{f}\left(\mathit{u}\right)\right)}{\mathit{\partial }\mathit{z}}=\mathit{p}.$$

(10)

Equation (10) is the PDE formulation which is solved using the CE/SE method. For the detailed derivation of the CE/SE method, the reader is referred to [12,17,18,19,20,21]. Here, only the final equations of the fully discretized scheme in space and time are presented. The state variable u at the j-th spatial point and n-th time instant is

$$\begin{gathered} {\mathit{u}}_j^n=\frac{1}{2}\left[{\mathit{u}}_{j-1/2}^{n-1/2}+{\mathit{u}}_{j+1/2}^{n-1/2}+{\mathit{s}}_{j-1/2}^{n-1/2}-{\mathit{s}}_{j+1/2}^{n-1/2}+\frac{∆t}{2}\left({\mathit{p}}_{j-1/2}^{n-1/2}+{\mathit{p}}_{j+1/2}^{n-1/2}\right)\right], \\ j=1, 2, \dots , N_z; n=1, 2, \dots , N_t, \end{gathered}$$

(11)

where ${\mathit{s}}_j^n$ is

$${\mathit{s}}_{\mathit{j}}^{\mathit{n}}=\frac{∆z}{\mathbf{4}}{\mathit{u}}_{\mathit{z},\mathit{j}}^n+\frac{∆\mathit{t}}{∆\mathit{z}}{\mathit{f}}_{\mathit{j}}^{\mathit{n}}+\frac{∆{\mathit{t}}^{\mathbf{2}}}{\mathbf{4}∆\mathit{z}}{\mathit{f}}_{\mathit{t},\mathit{j}}^{\mathit{n}}.$$

(12)

In the last equation, u_z and f_t are the discrete analogues of the derivatives $\partial \mathit{u}/\partial z$ and $\partial \mathit{f}/\partial t$. The current value of ${\mathit{u}}_j^n$ is calculated from the already available values of ${\mathit{u}}_{j\pm 1/2}^{n-1/2}, {\mathit{p}}_{j\pm 1/2}^{n-1/2}, {\mathit{u}}_{z,j\pm 1/2}^{n-1/2}, {\mathit{f}}_{j\pm 1/2}^{n-1/2},$ and ${\mathit{f}}_{t,j\pm 1/2}^{n-1/2}$ at the previous time instant. The different versions of the CE/SE method only differ in the way the spatial derivative u_z of the state variables is calculated. In the current work, the proposition of [21] is used.

$${\mathit{u}}_{\mathit{z},\mathit{j}}^{\mathit{n}}={\mathit{W}}_{\mathbf{0}}\left({\mathit{u}}_{\mathit{z}\mathbf{-},\mathit{j}}^{\mathit{n}},{\mathit{u}}_{\mathit{z}\mathbf{+},\mathit{j}}^{\mathit{n}}, \mathit{\alpha }\right).$$

(13)

The definition of the function W₀ is

$$\{\begin{array}{c}{\mathit{W}}_{\mathbf{0}}\left(\mathbf{0}, \mathbf{0},\mathit{\alpha }\right)=\mathbf{0}\\ {\mathit{W}}_{\mathbf{0}}\left({\mathit{x}}_{\mathbf{-}},{\mathit{x}}_{\mathbf{+}},\mathit{\alpha }\right)=\frac{{\left|{\mathit{x}}_{\mathbf{+}}\right|}^{\mathit{\alpha }}{\mathit{x}}_{\mathbf{-}}\mathbf{+}{\left|{\mathit{x}}_{\mathbf{-}}\right|}^{\mathit{\alpha }}{\mathit{x}}_{\mathbf{+}}}{{\left|{\mathit{x}}_{\mathbf{+}}\right|}^{\mathit{\alpha }}\mathbf{+}{\left|{\mathit{x}}_{\mathbf{-}}\right|}^{\mathit{\alpha }}}\end{array} , \left(\left|{\mathit{x}}_{\mathbf{+}}\right|\mathbf{+}\left|{\mathit{x}}_{\mathbf{-}}\right|>\mathbf{0}\right),$$

(14)

with $\mathit{\alpha }\ge \mathbf{0}$ [18]. In this study, $\mathit{\alpha }=\mathbf{1}$. ${\mathit{u}}_{\mathit{z}\mathbf{-},\mathit{j}}^{\mathit{n}}$ and ${\mathit{u}}_{\mathit{z}\mathbf{+},\mathit{j}}^n$ are the values of ${\mathit{u}}_{\mathit{z},\mathit{j}}^{\mathit{n}}$ calculated from left and right of the point with coordinates (j, n), respectively, and are given by

$${\mathit{u}}_{\mathit{z}\pm ,\mathit{j}}^{\mathit{n}}=\pm \frac{{\tilde{\mathit{u}}}_{\mathit{j}\pm \mathbf{1}/\mathbf{2}}^{\mathit{n}}-{\mathit{u}}_{\mathit{j}}^{\mathit{n}}}{∆\mathit{z}/\mathbf{2}},$$

(15)

where

$${\tilde{\mathit{u}}}_{\mathit{j}\pm \mathbf{1}/\mathbf{2}}^{\mathit{n}}={\mathit{u}}_{\mathit{j}\pm \mathbf{1}/\mathbf{2}}^{\mathit{n}-\mathbf{1}/\mathbf{2}}+\frac{∆\mathit{t}}{\mathbf{2}}{\mathit{u}}_{\mathit{t},\mathit{j}\pm \mathbf{1}/\mathbf{2}}^{\mathit{n}-\mathbf{1}/\mathbf{2}},$$

(16)

where u_t is the numerical analogue of $\partial \mathit{u}/\partial \mathit{t}$ and is

$${\mathit{u}}_{\mathit{t},\mathit{j}}^{\mathit{n}}=-{\mathit{f}}_{\mathit{z},\mathit{j}}^{\mathit{n}}+{\mathit{p}}_{\mathit{j}}^{\mathit{n}}.$$

(17)

Since the CE/SE method is a fully discrete explicit scheme, the Courant–Friedrichs–Lewy condition has to be fulfilled to guarantee numerical stability. This condition gives the relationship between the spatial and time step sizes through the Courant–Friedrichs–Lewy (CFL) number,

$$\mathrm{CFL}=\mathit{v}\frac{∆\mathit{t}}{∆\mathit{z}}\le \mathbf{1}.$$

(18)

In practice, we specify $∆z$ or $∆t$ and the CFL number and then calculate the other unspecified step $∆t$ or $∆z$.

4. Results

To compare the cell model (MOL) and CE/SE methods for the simulation of chromatographic processes in terms of computational efficiency, several examples were considered. All simulations were carried out on a desktop personal computer (PC) with Intel Core i7 6700 central processing unit (CPU) (3.4 GHz), with 8 GB random-access memory (RAM) on MathWorks MATLAB R2017a under the Microsoft Windows 8.1 operating system.

4.1. Single Column with the Ideal Equilibrium Model

First, the focus was placed on the equilibrium model without axial dispersion which can be used for the simulation of highly efficient columns. In practice, column efficiency is often specified in terms of number of theoretical stages N corresponding to the number of grid points of the first order finite volume MOL formulation. For N $\to \infty$ and a step input, discontinuities travel through the column during the loading of an empty bed for the linear and the Langmuir isotherms considered in this paper. For a linear isotherm, these are contact discontinuities, and, for the Langmuir isotherm, these are self-sharpening shock waves [2].

In practical situations, the number of theoretical plates can be high but will be finite. In the remainder, N is fixed and, therefore, so is the number of grid points of the MOL formulation. The number of grid points for the CE/SE method is adapted accordingly to meet the given column efficiency. For comparison, the analytical solution for N $\to \infty$ is also given.

The first scenario is concerned with the loading of a single column with a binary system described by Langmuir isotherms.

Table 1. Simulation parameters and reversed CE/SE method parameters for Example 1 (single-column binary process with Langmuir isotherms described by the ideal equilibrium model).

Quantity	Value	Quantity	Value	Quantity	Value
Lcol (m)	1.0	HA	2	HB	4
v (m·s−1)	0.1	bA (L·g−1)	2	bB (L·g−1)	4
ε	0.8	CFe,A (g·L−1)	0.9	CFe,B (g·L−1)	0.8
tsim (s)	10	Nt	501	$\overline{\mathrm{CFL}}$	0.4

summarizes all simulation parameters. Feed concentrations of the two components serve as boundary conditions for the PDE system Equation (6).

As already mentioned, the CE/SE method formulation is given by Equation (10), while the ideal chromatographic model is of type

$$\frac{\partial \left(\mathit{f}\left(\mathit{u}\right)\right)}{\partial t}+\frac{\partial \mathit{u}}{\partial z}=0.$$

(19)

For nonlinear isotherms, the model in Equation (19) can only be directly converted to the type in Equation (10) in special cases. A detailed discussion for Langmuir isotherms is given in Appendix A. A possible alternative is to interchange the spatial and the time coordinates, i.e., instead of propagating in time, the solution procedure propagates in space (Figure 1b). This particular use of the CE/SE method is called the reversed CE/SE method. The CFL number is then formulated as

$$\overline{\mathrm{CFL}}=\frac{1}{v}\frac{∆z}{∆t}.$$

(20)

The algorithm is as follows:

(i)

Specify the simulation time t_sim;

(ii)

Specify the number of time steps N_t, i.e., the number of conservation elements;

(iii)

Calculate the time step size $∆t$;

(iv)

Specify $\overline{\mathrm{CFL}}$ ($0<\overline{\mathrm{CFL}}\le 1)$;

(v)

Calculate the spatial step size $∆z$ from Equation (20);

(vi)

Continue with the reversed CE/SE method.

For the solution obtained via the MOL, the integration in time was done using the built-in MATLAB ODE solver ode45 which uses the Runge–Kutta method [25]. Doubled computational time was found using ode23, which is also a Runge–Kutta method but of lower order. As an alternative for stiff systems, the ode15s solver was also used but it was threefold slower than ode45. This indicates that the present system is nonstiff.

The concentration profiles along the column calculated by the two numerical methods are shown in Figure 2a for three time instants—0, 5, and 10 s. Additionally, the analytical solution of the system in Equation (4) obtained via the method of characteristics for N$\to \infty$ is shown [26]. It can be seen that both MOL and reversed CE/SE methods resolved the concentration shocks with the same number of points. For MOL, the number of points refers to the number of finite volumes in space, while, for the reversed CE/SE method, the number of points refers to the number of elements in time. The corresponding number of elements in space followed from the given number of elements in time and the given $\overline{\mathrm{CFL}}$ number from Equation (20). The computational time of the MOL (ode45) simulation was 7.85 s, which was the reference time for all of the simulations in the current work. The reversed CE/SE method was much more efficient (approximately 34-fold) as shown in Figure 2b, even with the chosen $\overline{\mathrm{CFL}}$ number of 0.4, which is rather conservative.

4.2. Single Column with the LDF Model

Unfortunately, the reversed formulation of the CE/SE method, as well as the other methods discussed in Appendix A, can only be used for single-column simulations because the feed concentrations, i.e., the boundary conditions are known a priori at all time steps. However, this is not true anymore for SMB processes and control problems. That is why, for systems with nonlinear Langmuir isotherms, the mass-transfer model in Equation (1) was used. The LDF chromatographic model corresponded to Equation (10), and the original formulation of the CE/SE method could be used for simulations. Interchanging of space and time coordinates was not required. It is worth noting that the LDF model represents a system of semi-linear transport equations compared to the equilibrium model, which is quasi-linear and, therefore, represents a nonlinear system of transport equations [27]. Thus, the LDF model has a simpler mathematical structure, but the number of equations is doubled.

The parameters for this example were the same as those for Example 1 presented in Table 1. Additional parameters were the mass transfer coefficients. Here, the same value was assumed for both components. The value of k_m and the CE/SE method parameters are given in

Table 2. Simulation parameters and CE/SE method parameters for Example 2 (single-column binary process with Langmuir isotherms described by the linear driving force (LDF) model).

Quantity	Value	Quantity	Value
km,A (s−1)	10	Nz	101
km,B (s−1)	10	CFL	0.4

. The values of the mass transfer coefficients were chosen to be high enough such that physical mass transfer was negligible.

The results of the CE/SE and MOL simulations are shown in Figure 3a together with the analytical solution for $N\to \infty$ and $k_{m,i}\to \infty$. It can be seen that the orginal CE/SE method formulation could resolve the concentration fronts very well with a small number of spatial points, i.e., 101 points, whereas MOL needed more than 2000 points to achieve a similar accuracy. In both cases, the number of points refers to the spatial discretization. Therefore, the computational times of the two methods differed significantly (Figure 3b), and the CE/SE method was more than 55 times faster than MOL. Again, ode45 was the solver used for MOL. For the LDF model, computational times for ode45 and ode23 were found to be very similar. In particular, ode45 or ode23 was about 280-fold faster than ode15s, indicating that the LDF model was considerably less stiff than the equilibrium one. Therefore, in the remainder of the study, unless stated otherwise, ode45 was used for MOL.

For lower values of the mass transfer coefficients, k_m,_i concentration fronts were less steep compared to those in Figure 3a. In consequence, a lower number of grid points was required to resolve the fronts in both cases (CE/SE and MOL), leading to lower computational times.

4.3. Binary SMB Process with the LDF Model

The dynamical model of the SMB processes includes the PDE system in Equation (1) to describe each chromatographic column, as well as the nodal material balances which define the boundary conditions for the system in Equation (1). These equations are shown below.

• − Desorbent node

$$Q_{De}=Q_I-Q_{IV};$$

(21)

$$Q_{\mathrm{De}}{C}_{i,\mathrm{De}}=Q_{\mathrm{I}}{C}_{i,\mathrm{I}}^{\mathrm{in}}-Q_{\mathrm{IV}}{C}_{i,\mathrm{IV}}^{\mathrm{out}}, i=1,\dots , N_{\mathrm{comp}}.$$

(22)

• − Extract node

$$Q_{Ex}=Q_I-Q_{II};$$

(23)

$$C_{i,\mathrm{Ex}}=C_{i,\mathrm{I}}^{\mathrm{out}}=C_{i,\mathrm{II}}^{\mathrm{in}}, i=1,\dots , N_{\mathrm{comp}}.$$

(24)

• − Feed node

$$Q_{Fe}=Q_{III}-Q_{II};$$

(25)

$$Q_{\mathrm{Fe}}{C}_{i,\mathrm{Fe}}=Q_{\mathrm{III}}{C}_{i,\mathrm{III}}^{\mathrm{in}}-Q_{\mathrm{II}}{C}_{i,\mathrm{II}}^{\mathrm{out}}, i=1,\dots , N_{\mathrm{comp}}.$$

(26)

• − Raffinate node

$$Q_{Ra}=Q_{III}-Q_{IV};$$

(27)

$$C_{i,\mathrm{Ra}}=C_{i,\mathrm{III}}^{\mathrm{out}}=C_{i,\mathrm{IV}}^{\mathrm{in}}, i=1,\dots , N_{\mathrm{comp}}.$$

(28)

The process configuration for the binary SMB process is shown in Figure 4a. In this case study, each zone had only one chromatographic column. The simulation parameters and the operating conditions are presented in

Table 3. Simulation parameters for Example 3 (binary SMB process with Langmuir isotherms described by the LDF model).

Quantity	Value	Quantity	Value	Quantity	Value	Quantity	Value
Lcol (m)	0.5	mI	5	HA	2	HB	4
Dcol (m)	0.02	mII	1.8	bA (L·g−1)	0.2	bB (L·g−1)	0.4
ε	0.8	mIII	2.8	km,A (s−1)	10	km,B (s−1)	10
Dpipe (m)	0.002	mIV	1.3	CFe,A (g·L−1)	0.9	CFe,B (g·L−1)	0.7
tsw (s)	40

. The operating point, i.e., the m values, was close to the optimal operating point for this system according to the triangle theory [4]. The m values are dimensionless flowrate ratios of liquid and solid phases. For the SMB process, these are defined as

$$m_p=\frac{\mathrm{Liquid}\mathrm{phase}\mathrm{flowrate}}{\mathrm{Solid}\mathrm{phase}\mathrm{flowrate}}=\frac{Q_{\mathrm{int},p}{t}_{sw}-\epsilon V_{\mathrm{col}}}{\left(1-\epsilon \right)V_{\mathrm{col}}},$$

(29)

where t_sw is the switching time (s), V_col is the column volume (m³), and p is the zone index. From these values, the internal flowrates Q_int (m³·s⁻¹) in each zone of the SMB plant could be calculated.

After several column switchings, a cyclic steady state (CSS) was reached. In this particular example, this cyclic steady state was reached after 25 cycles. Figure 4c,d show the concentration profiles of the two components inside the SMB plant during this startup. Dashed–dotted curves are the profiles at the beginning of every cycle just after the column switching, while the solid curves correspond to the profiles at end of the cycle. With the chosen values for the isotherm parameters H_i and b_i, the component A (blue) had lower affinity to the solid phase and flowed out of the raffinate port, while the component B (red) had higher affinity to the solid phase and flowed out from the extract port.

For the simulation with the CE/SE method, the value of the CFL number was selected as 0.4. To satisfy the CFL condition in each zone of the SMB, the value of the time step $∆t$ size was calculated from the highest liquid velocity in the system, i.e., in zone I.

$$∆t=\mathrm{CFL}\frac{∆z}{\max v_p}.$$

(30)

The number of spatial points was 101 per column, and the CSS was reached after 47% relative computational time. To achieve similar accuracy with the MOL, 10 times more spatial discretization points per column were needed. This led to the much slower performance of the MOL simulation (Figure 4b).

For MOL in Figure 4b,d, ode45 was used. Again, the performance of ode23 was similar (1002% compared to 1444.6% for ode45). In contrast, ode15s needed 580-fold more computational time than the simulation with ode45.

4.4. Ternary Center-Cut Eight-Zone SMB Process with Linear Isotherms and the Ideal Equilibrium Model

Ternary SMB processes, called center-cut separations, play an important role in biotechnology and pharmaceutical industries for the isolation of desired key components that have medium affinity to the solid phase in comparison to the other fractions in the mixture. There are several configurations for ternary center-cut separation—cascade of two SMB units, eight-zone SMB, Japan Organo process, etc. [28]. Eight-zone SMB processes can be used with raffinate or with extract recycle [23]. In the present work, an eight-zone SMB process with raffinate recycle and both linear and Langmuir isotherms was investigated.

Since, for systems with linear isotherms, the ideal equilibrium model Equation (4) can be converted directly to form Equation (10), it was used as a first example. The process configuration is shown in Figure 5a. Again, each zone had one column. The nodal material balances for the eight-zone SMB with raffinate recycle are as follows:

• − First desorbent node

$$Q_{De1}=Q_{I1}-Q_{IV2};$$

(31)

$$Q_{\mathrm{De}1}{C}_{i,\mathrm{De}1}=Q_{\mathrm{I}1}{C}_{i,\mathrm{I}1}^{\mathrm{in}}-Q_{\mathrm{IV}2}{C}_{i,\mathrm{IV}2}^{\mathrm{out}}, i=1,\dots , N_{\mathrm{comp}}.$$

(32)

• − First extract node

$$Q_{Ex1}=Q_{I1}-Q_{II1};$$

(33)

$$C_{i,\mathrm{Ex}1}=C_{i,\mathrm{I}1}^{\mathrm{out}}=C_{i,\mathrm{II}1}^{\mathrm{in}}, i=1,\dots , N_{\mathrm{comp}}.$$

(34)

• − First feed node

$$Q_{Fe1}=Q_{III1}-Q_{II1};$$

(35)

$$Q_{\mathrm{Fe}1}{C}_{i,\mathrm{Fe}1}=Q_{\mathrm{III}1}{C}_{i,\mathrm{III}1}^{\mathrm{in}}-Q_{\mathrm{II}1}{C}_{i,\mathrm{II}1}^{\mathrm{out}}, i=1,\dots , N_{\mathrm{comp}}.$$

(36)

• − First raffinate node

$$Q_{Ra1}=Q_{III1}-Q_{IV1};$$

(37)

$$C_{i,\mathrm{Ra}1}=C_{i,\mathrm{III}1}^{\mathrm{out}}=C_{i,\mathrm{IV}1}^{\mathrm{in}}, i=1,\dots , N_{\mathrm{comp}}.$$

(38)

• − Second desorbent node

$$Q_{De2}=Q_{I2}-Q_{IV1};$$

(39)

$$Q_{\mathrm{De}2}{C}_{i,\mathrm{De}2}=Q_{\mathrm{I}2}{C}_{i,\mathrm{I}2}^{\mathrm{in}}-Q_{\mathrm{IV}1}{C}_{i,\mathrm{IV}1}^{\mathrm{out}}, i=1,\dots , N_{\mathrm{comp}}.$$

(40)

• − Second extract node

$$Q_{Ex2}=Q_{I2}-Q_{II2};$$

(41)

$$C_{i,\mathrm{Ex}2}=C_{i,\mathrm{I}2}^{\mathrm{out}}=C_{i,\mathrm{II}2}^{\mathrm{in}}, i=1,\dots , N_{\mathrm{comp}}.$$

(42)

• − Second feed node

$$Q_{Fe2}=Q_{Ra1};$$

(43)

$$C_{i,\mathrm{Fe}2}=C_{i,\mathrm{Ra}1}, i=1,\dots , N_{\mathrm{comp}}.$$

(44)

• − Second raffinate node

$$Q_{Ra2}=Q_{III2}-Q_{IV2};$$

(45)

$$C_{i,\mathrm{Ra}2}=C_{i,\mathrm{III}2}^{\mathrm{out}}=C_{i,\mathrm{IV}2}^{\mathrm{in}}, i=1,\dots , N_{\mathrm{comp}}.$$

(46)

Table 4. Simulation parameters for Example 4 (ternary center-cut eight-zone SMB process with linear isotherms described by the ideal equilibrium model).

Quantity	Value	Quantity	Value	Quantity	Value	Quantity	Value	Quantity	Value
Lcol (m)	0.5	mI,1	2.55	mI,2	1.82	HA	1.1	CFe,A (g·L−1)	0.9
Dcol (m)	0.02	mII,1	1.57	mII,2	1.22	HB	1.7	CFe,B (g·L−1)	0.8
ε	0.75	mIII,1	2.19	mIII,2	2.55	HC	2.5	CFe,C (g·L−1)	0.7
Dpipe (m)	0.002	mIV,1	0.86	mIV,2	1.01	tsw (s)	60

summarizes the simulation parameters and the operating conditions, which were taken from [22]. Component C (red) had the highest affinity to the solid phase and flowed out through the first extract port together with some amount of the middle component B (green). The outlet mixture from the first raffinate port, which contained the components A (blue, the component with the lowest affinity to the solid phase) and B, was fed to the second subunit where it was separated. Component B flowed out from the second extract port, while component A flowed out from the second raffinate port. Simulation results show that CSS was reached after 32 cycles, and the concentration profiles are shown in Figure 5c,d. For the CE/SE method simulation, the CFL number was 0.8 and the number of spatial discretization points per column was 101. For the MOL, 401 such points were needed per column, and the mass matrix for systems with linear isotherms is constant and, therefore, was calculated once at the beginning of the simulation, leading to a very fast computational time. As result, the CE/SE method was only 11-fold faster than MOL (Figure 5b).

Finally, it is worth mentioning that, for linear isotherms and for N $\to \infty$, a new analytical solution approach is available, which was applied to binary and ternary center-cut SMB processes and which is even orders of magnitude faster than the CE/SE method [29].

4.5. Ternary Center-Cut Eight-Zone SMB Process with Langmuir Isotherms and the LDF Model

Lastly, a study of a ternary system with Langmuir isotherms was performed. The process configuration was the same as in Figure 5a, i.e., ternary eight-zone center-cut SMB with raffinate recycle. The nodal material balances were those in Equations (31)–(46). Physical parameters of the chromatographic columns are presented in Table 4 (first column), while the operating conditions and the parameters of the participating components are listed in

Table 5. Simulation parameters for Example 5 (ternary center-cut eight-zone SMB process with Langmuir isotherms described by the LDF model).

Quantity	Value	Quantity	Value	Quantity	Value	Quantity	Value	Quantity	Value
mI,1	2.55	mI,2	2.10	HA	1	HB	2	HC	2.5
mII,1	1.893	mII,2	0.928	bA (L·g−1)	1	bB (L·g−1)	2	bC (L·g−1)	2.5
mIII,1	1.90	mIII,2	1.99	km,A (s−1)	10	km,B (s−1)	10	km,C (s−1)	10
mIV,1	0.915	mIV,2	0.85	CFe,A (g·L−1)	0.5	CFe,B (g·L−1)	5.0	CFe,C (g·L−1)	1.5
tsw (s)	40

. The operating conditions were calculated with the methodology presented in [30] for a true moving bed process (TMB) with little adjustments to take the SMB process into account; to the best of our knowledge, this is the first simulation of ternary eight-zone SMB with Langmuir isotherms in the nonlinear concentration range. Calculations for the nonlinear case were based on the LDF model.

Concentration profiles along the SMB plant are shown in Figure 6a,b. Target component B was obtained as a pure component from the extract port of the second SMB subunit. The cyclic steady state was reached after around 60 cycles. The CFL number for the CE/SE method was 0.4, and the number of spatial discretization points was 101 per column. On the other hand, MOL needed 10-fold more spatial points per column to achieve a similar accuracy to CE/SE. This led to enormous computational time for the MOL method (around 2.5 h) in comparison to the CE/SE simulation (Figure 6c), which is 350-fold greater than that by CE/SE.

5. Discussion and Conclusions

This paper focused on the numerical simulation of chromatographic processes using the explicit space–time CE/SE method. For a systematic evaluation of the method, different models with linear and nonlinear isotherms and different process configurations including single-column, binary, and ternary eight-zone SMB processes were investigated.

For the first time, the application of the CE/SE method to the popular equilibrium model was also considered. It was shown that application of the CE/SE method requires reformulation of the equilibrium model equations. For linear isotherms, this reformulation is easy to achieve and straightforward. For nonlinear isotherms, reformulation is more involved. In particular, two different approaches were proposed, namely, (i) the reversed space–time formulation, which can be applied to any nonlinear adsorption isotherm, and (ii) an inversion-based approach, which applies only to the important but limited class of nonlinear Langmuir isotherms. In either case, reformulation is only possible for single-column processes. For the reversed method, the solution procedure propagates in space instead of time. Therefore, the temporal evolution of the boundary conditions needs to be known before the start of the simulation, which is not possible for multicolumn processes with recycle. For the inversion-based method, discussed in detail in Appendix A, an adjusted time was introduced, which depended on the linear velocity of the column and would be different from column to column in a multicolumn process, thus preventing simultaneous integration of multicolumn processes in time. Consequently, the LDF model was used for the simulation of nonlinear multicolumn processes since reformulation is not required for the LDF model.

In all cases, the CE/SE method was shown to be much more efficient than the popular cell model, which represents a method of lines approach using a first-order finite volume discretization (MOL) in combination with established time integrators from MATLAB such as the Runge–Kutta method ode45. Computational efficiency of the two methods was measured in terms of computational times for same resolution of the concentration profiles. The largest reduction in computation times was found for processes with steep fronts and with nonlinear Langmuir adsorption isotherms. In particular, the CE/SE method was found to be about 350-fold faster than MOL with ode45 for a ternary eight-zone SMB process with raffinate recycle for a challenging center-cut separation. This opens new possibilities for an efficient computational evaluation of different processes for ternary center-cut separations with nonlinear adsorption isotherms. The main focus in this field has so far been on processes with linear adsorption isotherms described by the true moving bed approximation [23].

In our future work, we also intend to use the CE/SE method for model-based control of binary and ternary SMB processes following the approach in [31,32,33].