E ﬃ cient Simulation of Chromatographic Processes Using the Conservation Element / Solution Element Method

: Chromatographic separation processes need e ﬃ cient 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 di ﬀ erential equations. Numerical stability of this method is guaranteed when the Courant–Friedrichs–Lewy condition is satisﬁed. To investigate the accuracy and e ﬃ ciency of this method, it is compared with the classical cell model, which corresponds to a ﬁrst-order ﬁnite volume discretization using a method of lines approach (MOL). The evaluation is done for di ﬀ erent models, including the ideal equilibrium model and a mass transfer model for di ﬀ erent adsorption isotherms—including linear and nonlinear Langmuir isotherms—and for di ﬀ erent 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.


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.

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.
where ε is the column void fraction, v is the liquid phase velocity (m·s −1 ), D ax is the axial dispersion coefficient (m 2 ·s −1 ), k m is the mass transfer coefficient between the two phases (s −1 ), C is the liquid-phase concentration (g·L −1 ), q is the solid-phase concentration (g·L −1 ), q* is the solid-phase concentration at the interphase boundary in equilibrium with the liquid phase (g·L −1 ), 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 → 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: 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 −1 ), and r is the component index. In the limiting case where the mass transfer is instantaneous (i.e., k m,i → ∞ and q i, k → q * i, k ) and there is negligible dispersion, the ideal equilibrium model [1] is obtained.
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.
Boundary conditions (BCs) are of the Dirichlet type, 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.

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 Dax → 0 can be presented in the following vector form (the indices are omitted for simplicity): 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: 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 where is In the last equation, uz and ft are the discrete analogues of the derivatives ⁄ and ⁄ .
The current value of is calculated from the already available values of ± / / , ± / / , , ± / / , ± / / , and , ± / / at the previous time instant. The different versions of the The model in Equation (1) with D ax → 0 can be presented in the following vector form (the indices are omitted for simplicity): 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: 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 where s n j is In the last equation, u z and f t are the discrete analogues of the derivatives ∂u/∂z and ∂f/∂t. The current value of u n j is calculated from the already available values of u n−1/2 j±1/2 , p n−1/2 j±1/2 , u n−1/2 z,j±1/2 , f n−1/2 j±1/2 , and f n−1/2 t,j±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. u n z,j = W 0 u n z−,j , u n z+,j , α .
The definition of the function W 0 is with α ≥ 0 [18]. In this study, α = 1. u n z−,j and u n z+,j are the values of u n z,j calculated from left and right of the point with coordinates (j, n), respectively, and are given by where u n j±1/2 = u n−1/2 j±1/2 + ∆t 2 u n−1/2 t,j±1/2 , where u t is the numerical analogue of ∂u/∂t and is 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, In practice, we specify ∆z or ∆t and the CFL number and then calculate the other unspecified step ∆t or ∆z.

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.

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 → ∞ 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 → ∞ 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 summarizes all simulation parameters. Feed concentrations of the two components serve as boundary conditions for the PDE system Equation (6). 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 As already mentioned, the CE/SE method formulation is given by Equation (10), while the ideal chromatographic model is of type 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 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 CFL (0 < CFL ≤ 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 → ∞ 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 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 CFL number of 0.4, which is rather conservative.

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 km and the CE/SE method parameters are given in Table 2. The values of the mass transfer coefficients were chosen to be high enough such that physical mass transfer was negligible. 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 The results of the CE/SE and MOL simulations are shown in Figure 3a together with the analytical solution for N → ∞ and km,i → ∞. 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

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. The values of the mass transfer coefficients were chosen to be high enough such that physical mass transfer was negligible. 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 The results of the CE/SE and MOL simulations are shown in Figure 3a together with the analytical solution for N → ∞ and k m,i → ∞ . 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. 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, km,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.

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.
− Feed node   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.

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 − Extract node − Feed node − Raffinate node 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. 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 where t sw is the switching time (s), V col is the column volume (m 3 ), and p is the zone index. From these values, the internal flowrates Q int (m 3 ·s −1 ) in each zone of the SMB plant could be calculated.

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 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 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.
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.

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  Finally, it is worth mentioning that, for linear isotherms and for N → ∞, 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]. (33) − First raffinate node − Second desorbent node − Second feed node Table 4 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). 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 Finally, it is worth mentioning that, for linear isotherms and for N → ∞ , 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].

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. 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. 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 m I, 1 2.55 m I,2 2.10 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.

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. 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.

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

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].  Acknowledgments: Valentin Plamenov Chernev would like to thank to Heorhii Marhiiev from Donetsk National Technical University in Pokrovsk (DonNTU), Ukraine for the fruitful discussions related to the CE/SE method during his stay in Magdeburg.

Conflicts of Interest:
The authors declare no conflict of interest. The funders had no role in the design of the study; in the collection, analysis, or interpretation of data; in the writing of the manuscript, or in the decision to publish the results.

Appendix A. Direct Conversion of the Equilibrium Model to the Form Given by Equation (10)
As already mentioned above, the CE/SE method applies to PDEs of the type shown in Equation (10). In contrast to this, the equilibrium model of chromatography is given by with u = C and f(u) = εC + (1 − ε)q(C). In the case of a linear isotherm, shown in Equation (2), its derivative is constant, and Equation (A1) can be easily rearranged to the form in Equation (10). For nonlinear isotherms, an inversion of the function f(u) is required and the corresponding system of PDEs can then be solved for f instead of u according to ∂f ∂t which is again equivalent to the form given by Equation (10). For the single-component Langmuir isotherm, inversion leads to the solution of a single quadratic equation, which can be solved explicitly for the unique positive solution to bring the model to the form given by Equation (10). For the multicomponent Langmuir isotherm, inversion leads to the solution of a system of quadratic equations, which is much more involved. However, the problem can be simplified by transformation of the equilibrium model to a simpler form by using an adjusted time τ = vt − z [34]. With this adjusted time, the equilibrium model reads ∂(q(C)) ∂τ In a second step, the well-known inverse of the Langmuir isotherm [2] (p. 252) is applied, and Equation (A3), together with Equation (A4), is solved for q. For the Langmuir isotherm, this is always possible for any number of components.