Development of a fourth-order compact finite difference scheme for simulation of simulated-moving-bed process

A fourth-order compact finite difference scheme was developed to solve the model equation of simulated moving bed, which has a boundary condition that is updated along the calculation process and cannot be described as an explicit function of time. Two different methods, direct method and pseudo grid point method, were proposed to deal with the boundary condition. The high accuracy of the two methods was confirmed by a case study of solving an advection-diffusion equation with exact solution. The developed compact finite difference scheme was then used to simulate the SMB processes for glucose-fructose separation and enantioseparation of 1,1′-bi-2-naphtol. It was found that the simulated results fit well with the experimental data. Furthermore, the developed method was further combined with the continuous prediction method to shorten the computational time and the results showed that, the computational time can be saved about 45%.

The simulated moving bed (SMB) is a continuous preparative chromatography technique that has been widely used in various industries including petroleum 1,2 , food 3,4 , pharmaceutical 5-10 and biotechnology [11][12][13] to separate structurally similar compounds. To obtain an optimal operation conditions, a mathematical model is usually desired as it plays an essential role in the design and optimization of SMB process. Among them, the transport-dispersive model with linear driving force (LDF) approximation for describing mass transfer resistance between the mobile phase and the solid phase is believed to be the most widely used one with the following model equation 14,15 : e This model equation is solved for every switching period to obtain the concentration profiles in the liquid and solid phase at any time. The initial condition of each switching period depends on the concentration profile in the columns at the end of the former switching period, and for the first switching period, the initial condition is: The boundary conditions at the column inlet (x = 0) and outlet (x = L) are: The mathematical SMB model equations have been solved using various methods such as the finite element 16 , finite difference 17 , finite volume (especially WENO) 14,[18][19][20][21][22] , wavelet collocation 15 , as well as the space time conservation element and solution element (CE/SE) method 14,23 . However, no matter the method used, the accuracy and efficiency are usually in a contradict manner, i.e. a method with higher accuracy generally requires extensive calculations and vice verse. Therefore, developing a new method with high accuracy and efficiency, or at least making a suitable compromise between them is still a major issue for the design and optimization of the SMB processes.
The compact finite difference scheme (CFDS) attracts a great attention in the last twenty years [24][25][26] . In conventional finite difference scheme, the central difference with 2nd order precision needs 3 grid points. To construct a 4th order approximation, at least 5 grid points are necessary 27 , which would inevitably increase computation time greatly and add complexity for handling the boundary condition. However, CFDS can give a 4th order scheme by using only 3 grid points, which is achieved by coupling with the original governing equation 28 . That means CFDS combines a higher precision and fewer computation. However, the situation that derivation of CFDS needs governing equation makes it difficult to generalize. Usually the CFDS is constructed in a case by case manner, which is prohibitive for containing tedious mathematical treatment. Nevertheless, CFDS has been successfully used to solve numerous equations including convection-diffusion equation [29][30][31][32][33][34][35] , heat equation 25,36,37 , Gross-Pitaevskii equation 38 and Helmholtz equation 39 29,32 . But it should be addressed that in their methods, the boundary condition is an explicit function of time, and the derivative of this function is required to deal with the boundary condition. While for the SMB equations, as shown in Eqs. 1-5, the value of c in in Eq. 4 is updated continuously with the calculation and as a consequence, the function of c in verse t cannot be explicitly described, which makes it impossible to obtain the derivative of c in . So the methods reported in literature 29,32 cannot be used to solve the SMB model equations.
To develop a 4th-order CFDS that can be used to solve the model equation of SMB, in the present work we compared two different methods for the handling of the boundary conditions. The accuracy of the developed method was examined by solving an equation with exact solution, and the computational efficiency was compared with the central difference scheme as well as the space-time solution element/conservation element method.

The Fourth-Order Compact Finite Difference Scheme
By discretization on spatial domain, the Eqs. 1-2 can be transformed to a set of ordinary differential equations with time as an independent variable that can be solved by traditional methods, such as the "odeint" method in Scipy modula of Python. We thus first present the compact finite difference scheme for the spatial derivatives by considering the partial differential equation as follows: The spatial domain [0, L] was discretized as follows: i with a constant step size h = L/(n − 1). For the interior points, Mohebbi and Dehghan 32 have derived the 4th order CFDS as below (the derivation can also be found in the Supplemental Materials): To construct the 4th order equations for the boundary points x 1 and x n , the forward and backward difference schemes were used by denoting as the forward and backward difference schemes of order j for the ith derivatives of c(x, t) about x, respectively. The formulas for i = 1, 2 and j = 1, 2 can be easily obtained from the Taylor expansions of c(x, t) at the points adjacent to x i and are listed in Table 1, some of them can also be obtained from literature 27 .
Two methods for handling of the boundary points were proposed here. In the first method, the derivatives at the boundary points were approximated by the forward or backward finite difference schemes directly thus named as "direct method". While in the second method, a pseudo grid point (x 0 ) was assumed to be exist at the left side of the boundary point x 1 , and on the contrary, a pseudo point (x n+1 ) is located at the right side of x n . So the second method was named as "pseudo grid point method". Direct method. The equation dealing with the left boundary was deduced first. With the Taylor expansions of c(x 2 , t) and c(x 3 , t), the first derivative of c(x, t) at x 1 can be expressed as:   Using Eqs. 12-13, and substituting the derivatives by the forward difference schemes (Table 1) with appropriate order, the following equations are obtained: Combining Eqs. 14-15, Eq. 11 and the boundary condition in Eq. 7 gives: x a x a x a x a Similarly, for the right boundary conditions, it can be obtained: x a x a x n a x a Pseudo grid point method. For the left boundary point, it is assumed that there exists a pseudo grid point x 0 at the left hand of x 1 . Then x 1 can be looked as an interior point, and by using Eq. (10) the following equation was obtained:

Symbol Derivative Truncation error Formula
∆ c x t ( , ) www.nature.com/scientificreports www.nature.com/scientificreports/ It should be emphasized that, in the right hand of Eq. (18), only the forward difference schemes were used for approximation of the derivatives of f(x 1 , t). This is because the partial difference of c about time was included in f, ∂ ∂ c x t t ( , )/ is meaningless at x 0 . Nevertheless, the replacement of central difference schemes (δ x 1 and δ x 2 ) by forward difference schemes (∆ x 2 1 and ∆ x 2 2 ) will not lead to the loss of accuracy, because all these four difference schemes have an accuracy of 2nd order.
In Eq. (18), c(x 0 , t) is implicitly included in the operators δ x 1 and δ x 2 . It needs be eliminated by using the boundary condition Eq. (7). For this purpose, the first order and third order of c(x 1 , t) are expressed as follows: where Eq. (19) can be obtained by using the Taylor expansions at c(x 1 , t) and c(x 2 , t), while Eq. (20) is similar to Eq. (14) except that the operator δ x 2 is used to substitute ∆ x 2 2 , both of them have a truncation error of O(h 2 ). Combining Eqs. 19-20, Eq. 11 and omitting the error term gives: x a x a By combining Eq. 18 and Eq. 21 to eliminate c x t ( , ) 0 , one obtains: 2 ( , ) 12 The right boundary conditions can be handled in a similar way, except that the backward difference schemes are used, the finally obtained equation is as follows: is as follows: And the boundary conditions are: While using the "direct method" for handling the boundary conditions, Eqs. 24-26 are introduced into Eq. 11 and Eqs. [16][17]. Then the SMB model equations are transformed to a system of ordinary equations as follows:   If the pseudo grid point method was used to treat the boundary conditions, the resulted equations have the same pattern as Eqs. 27-28, but the vectors g(t), A and B are different from Eqs. 32-34, which are as follows:  The resulted ordinary differential equations are solved by the method of "odeint" in Scipy modula of python 3.6.4. All the calculations in this work were conducted on a PC with CPU 3.60 GHz, RAM 16.0 GB.

Advection-diffusion equation with analytical solution.
To test the accuracy of the CFDS developed in this work, the transient one-dimensional advection-diffusion equation was first solved as follows 29,32,40 : with initial condition: For this equation, the ordinary equation Eq. 27 is simplified to: in which, while using the "direct method" to treat the boundary conditions,  when using the pseudo grid point method for treating the boundary conditions, A is the same as in Eq. 36, and B is as follows: The advection-diffusion equation was solved by the CFDS with two different methods for handling the boundary conditions. It was found that both methods have a high accuracy. As an example, the solutions obtained by the direct method at final time t = 1 with v = 0.3, h = 1/2 8 , and different Pe numbers (Pe = vL/D a ) are shown in Fig. 1. It can be observed that the calculated concentration profiles are in a good accordance with the exact solutions, even for a higher Peclet number of 1000 (Fig. 1c). When using the pseudo grid point method to handle the boundary conditions, the solution has no significant difference compared with that shown in Fig. 1. Then the maximum error was evaluated with the maximum difference between the concentrations obtained from CFDS and analytical solution while x ranging from 0 to 1, it was found appear at x = 0.5. The maximum errors of the two methods were shown in Fig. 2. Of course, for a fixed h value, the maximum error will increase with the increase of Pe number. This is reasonable because a higher Pe number results in a sharper peak, and as a consequence, a denser mesh (smaller h value) along with the x-axis is usually needed. Meanwhile, the two methods for handling the boundary condition have a similar accuracy at higher Pe numbers ( > 100), but at lower Pe number (say Pe = 10) the pseudo grid point method gives a higher accuracy than the direct method. So, the pseudo grid point method is recommended and hereinafter was used for further analysis. However, it should be noted that the CFDS presented in this work has a very high accuracy at lower Pe numbers, such as the maximum errors at Pe = 10 shown in Fig. 2 were 2.4×10 −7 and 2.6×10 −8 for the direct and pseudo grid point methods, respectively. So even the direct method is accurate enough for most of the applications.
To examine the convergence rate of the method developed above, the advection diffusion equation was solved at h = 1/2 6 , 1/2 7 , 1/2 8 , and 1/2 9 . For a method with n-order accuracy, the solution error (E) at any point is proportional to h n and so the date log(E) vs. log(h) should be asymptotic to a straight line with slope n. Therefore, the solution errors at different h are shown in a double logarithmic plot (Fig. 3) by measuring the error at x = 0.5, which is the maximum error along the x-axis. It can be seen that the points have a good linear relationship. By linear regression, the slope was found to be 3.996, which thus confirms that the method developed in this work has a convergence rate of 4th-order. As a comparison, the central differential scheme 41 was also used to solve this advection diffusion equation at different h values, and as expected, a convergence rate of 2-order was confirmed with a slope of 2.012 as shown in Fig. 3.
We also used the CFDS that was developed by Cao et al. 29 to solve the equations defined in Eqs. 38-41. The results showed that the maximum errors at different h values are the same as that of our method (Fig. 3) thus suggesting that the two methods have a similar accuracy. Because in Cao′s method, the vector A in Eq. 33 is tri-diagonal (for detailed information, please refer to the original paper), this makes their methods favorable for Scientific RepoRtS | (2020) 10:7820 | https://doi.org/10.1038/s41598-020-64562-8 www.nature.com/scientificreports www.nature.com/scientificreports/ the stability analysis. However, the differential function of boundary conditions, i.e. g′(t), is needed, thus makes Cao's method not possible for some applications, such as the SMB equation, where c in is included in the boundary conditions (Eq. 4). Because c in can only be obtained with calculation, its differential function and so g′(t) cannot be obtained. This is our aim in the present work to develop a compact finite difference scheme to solve SMB equation. Of course we believe our method can also be used for other applications because the differential function of boundary conditions is not needed. However, we have to point out that the matrix A in Eq. 33 is not tri-diagonal thus the stability analysis becomes difficult. Nonetheless, the stability and small stencil of CFDS have been widely www.nature.com/scientificreports www.nature.com/scientificreports/ accepted. For the detailed discussion of the stability while using CFDS for solving the convection-diffusion equations, please refer to the literature. 29,32,42,43 Numerical solution of simulated-moving-bed model equation. We then sought to use the compact finite difference scheme to solve the model equation of simulated moving bed. Two SMB processes, glucose-fructose separation and enantioseparation of 1,1'-bi-2-naphtol, were used as the case studies. The parameters of the two systems were obtained from literature 4,7,44,45 as summarized in Table 2. The adsorption isotherms of sugar separation process is nearly linear: = .  www.nature.com/scientificreports www.nature.com/scientificreports/ In the simulation, the step sizes are set to be h = L/64 and L/40 for the sugar separation and enantioseparation processes respectively 14,15,23 . Furthermore, the criterion for achieving the cyclic steady state (CSS) is set to be that the maximum difference between the concentrations in two consecutive iterations is lower than 10 −4 of feed concentration. It is found that the sugar separation process needs 91 switches for achieving CSS upon this criterion, while it is generally acknowledged that 80 switches are sufficient to reach CSS in this SMB process 4,15 . So it is believed that this criterion is appropriate and is used in the simulation of the two case studies.
The simulated concentration profiles at half of a switching period after reaching CSS are shown in Fig. 4. It can be seen that the simulation results fit the experimental data well for both the two SMB systems. For comparison, the space-time conservation element and solution element (CE/SE) method 23 was also used to simulate the  Table 2. Parameters of the simulated moving bed processes for sugar (fructose-glucose) separation and enantioseparation of 1,1′-bi-2-naphtol. www.nature.com/scientificreports www.nature.com/scientificreports/ two SMB processes. The results are listed in Table 3. Obviously, the products purities and recoveries obtained by different methods are similar. But the calculation time of CFDS is much longer than that of CE/SE method. This is reasonable because a higher accuracy is usually accompanied by an expensive computation. So, for the applications where the simulation efficiency is desired, the CE/SE method is preferred, while for the applications where high accuracy is desired, the CFDS method is a very good option.
In our previous work, a continuous prediction method was developed to improve the simulation efficiency of SMB process 46 . The key point of this continuous prediction method is the construction or prediction of the concentration at CSS using the concentrations in the last two iterations. The formula is as follows: Here, state variable u contains the concentrations of two solutes in the mobile phase and in the stationary phase, i.e. c A , c B , q A , and q B . The superscript m means the mth iteration and the acceleration factor θ is an empirical parameter. Then, the predicted state variable, u m CSS ( ) , is used as the initial value for the next iteration. Through this method, the iterations needed to achieve CSS can be reduced and thus the calculation time can be shortened. So, the continuous prediction method was used to improve the efficiency of CFDS in this work. The iteration numbers needed to reach CSS at different θ values are shown in Fig. 5. It can be seen that the computational efficiencies of the two SMB processes were both significantly improved with the aid of continuous prediction method. For the sugar separation process, an acceleration factor θ = 0.8 will lead to an iteration number of 50 for reaching CSS. While without the continuous prediction method used, i.e. θ = 0, the iteration number is 91. Accordingly, the calculation time was saved about 45% due to the CPU time is proportional to the iteration number 46 . As for the enantioseparation process, the calculation time can also be saved about 44% with θ = 0.7 (46 iterations) compared to θ = 0 (82 iterations). But to our surprise, the optimum θ values are much higher than the value (θ = 0.5) that we have previously recommended 23,46 . For comparison, the dependence of iteration number on θ when using CE/SE method was also shown in Fig. 5. The optimum θ values are 0.5 and 0.6 for the enantioseparation and sugar separation processes respectively, which is in consistence with our previous results 46 . The optimum θ values in CFDS are shifted 0.2 higher than that in CE/SE method for both the two SMB processes (0.5 to 0.7 for enantioseparation and 0.6 to 0.8 for sugar separation processes). The reason for this shift is unclear, but it is reasonable to draw a conclusion that the optimum θ value depends on the methods that are used to solve the model equation, even for the same SMB process.  Table 3. Simulation results and calculation times of CFDS and CE/SE methods for two SMB processes of sugar (glucose-fructose) separation and enantioseparation of 1,1′-bi-2-naphtol.

Conclusion
In the present work, a fourth-order compact finite difference scheme was successfully developed to solve the advection diffusion equations with Neumann boundary conditions, which do not need the boundary conditions to be differentiable function. Two different methods, direct method and pseudo grid point method, were proposed and used to handle the boundary conditions. The higher accuracy of the compact finite difference scheme was confirmed by a case study with analytical solution. It was found that the pseudo grid point method results in a higher accuracy than the direct method when the Pe number is low (such as 10). But for a moderate and high Pe numbers (such as 100 and 1000), the two methods give the same accuracy.
It should be pointed out that although the CFDS method can be used to solve the SMB model equations, the calculation time is much longer than the space-time conservation element and solution element method. This problem, however, can be solved by use of the continuous prediction method, which improves the calculation efficiency of CFDS significantly meanwhile saves the calculation time about 45%. In summary, it is tempting to speculate that the combination of 4th-order compact finite difference scheme described in the present work and the continuous prediction method can be used for most of the SMB processes and has a wide application potential, especially when a higher accuracy is desired.