Abstract

Translate

The efficiency of controlling the simulated moving bed (SMB) has long been a critical issue in the chemical engineering industry. Most existing research relies on finite element methods, which often result in lower control efficiency and are unable to achieve online control. To enhance control over the SMB process, this paper employs the Crank–Nicolson method to develop a discrete dynamical model. This approach allows for the investigation of system stability and convergence, fundamentally addressing the sources of error. During the discretization of partial differential equations (PDEs), two main types of errors arise: intrinsic errors from the method itself and truncation errors due to derivative approximations and the discretization process. Research indicates that for the former, the iterative process remains convergent as long as the time and spatial steps are sufficiently small. Regarding truncation errors, studies have demonstrated that they exhibit second-order behavior relative to time and spatial steps. The theoretical validation shows that the iteration works effectively, and simulations confirm that the finite difference method is stable and performs well with varying SMB system parameters and controller processes. This provides a solid theoretical foundation for practical, real-time online control.

Full text

Turn on search term navigation

Translate

1. Introduction

Simulated moving bed (SMB) technology is a sophisticated process employed for the separation and purification of chemical materials or mixtures. It operates on the principle of a moving bed, leveraging the dynamic flow characteristics of materials within the bed to achieve separation. Typically, SMB technology comprises a fixed bed alongside a system that mimics the behavior of a moving bed. By simulating the fluid flow and material movement within the bed, this technology replicates the separation processes of an actual moving bed. The essence of SMB technology lies in its ability to separate various components by modulating the direction or speed of fluid flow. This technology finds extensive application in industries such as chemical engineering, petrochemicals, and pharmaceuticals [1,2,3]. The process of the SMB is illustrated in Figure 1. The advent of SMB technology has significantly improved separation efficiency within the chromatography industry. Nonetheless, due to the numerous parameters inherent in industrial processes, optimizing these parameters through empirical experimentation is both time-consuming and costly. Technicians often rely on their experience to set initial control parameters for the SMB, a practice that frequently results in suboptimal operation [4,5,6].

Researchers have been focusing on analyzing SMB processes using mathematical models. The main models used are the general rate model and the equilibrium diffusion model. The general rate model is very detailed and closely matches real processes, but it is complex. On the other hand, the equilibrium diffusion model is simpler and works well when the concentrations of substances being separated are low [7]. This model can be further divided into different types based on the isotherm used: linear, Langmuir, modified Langmuir, and competitive Langmuir isotherms [8].

Irrespective of the model employed, these models typically comprise a set of partial differential equations for which analytical solutions are challenging to derive. The ultimate objective of SMB technology is to achieve effective control, necessitating a focus on the efficiency of the solution process for these models. To this end, it is crucial to investigate the convergence and stability of the computational models for the SMB.

Dunnebier et al. proposed and compared three different modeling and simulation methods, evaluating the solutions in terms of complexity, accuracy, and computational time. The newly developed closed-form solution demonstrated high accuracy and low computational time [9]. Andrade Neto et al. presented a self-adjusting nonlinear MPC method for the enantiomeric separation of Praziquantel in an analog moving bed apparatus that was based on the finite element method [10]. Muhammed et al. developed a rigorous dynamic model to compare the SMB with traditional distillation as the next best choice [11]. The SMB was optimized using a genetic algorithm, maximizing the Research Octane Number (RON) and gross margin as the objective functions. Z. Yan et al. applied the finite difference method to establish a matrix subspace iteration method, conducting numerical simulations on parameters such as isotherm, mass transfer resistance, cycle switching time, number of series-connected columns, column type, length, diameter, volume, feed concentration, and circulation flow rate to build a column movement model. This work laid the foundation for the optimal control of the column [12]. Leao et al. proposed a numerical solution process for simulating the transient and steady-state behavior of the simulated moving bed (SMB) based on a mathematical model [13]. Majeed et al. investigated the accuracy and efficiency of using orthogonal collocation methods to solve partial differential equations in separation engineering, with the numerical solutions provided being almost as accurate as those obtained from actual machinery [14]. Y. Kim et al. performed computational fluid dynamics (CFD) simulations using the finite element method to obtain concentration data and conducted experimental analysis using residence time distribution to support their findings [15]. The model integrates diffusion mass transfer mechanisms and dual-substrate Monod kinetics. WS-Lee et al. employed a mathematical dynamic model with limited experimental parameters, alongside a data-driven machine learning approach using real industrial data, to evaluate the performance of the SMB [16].

In the field of controller research, SMB systems exhibit a diverse array of approaches, with various controllers being designed specifically for different materials to be separated. Wei et al. proposed a control method based on a piecewise affine model. The simulation results demonstrate the feasibility of this approach, as it can fit the actual yield–production curve of the target compound and impurities to achieve the desired target yield [17]. Hoon et al. applied a data-driven deep Q-network, a model-free reinforcement learning method, to train near-optimal control strategies for the SMB process [18]. Natarajan and Lee proposed applying Repetitive Model Predictive Control (RMPC) to the SMB process, specifically for the chromatographic separation of phenylalanine–tryptophan mixtures [19]. Klatt et al. introduced a model-based optimization control scheme for SMB chromatographic separation and its application in the separation of fructose and glucose [20]. Carlos and Alain introduced a novel method for controlling the SMB chromatographic separation process [21]. Marrocos et al. proposed a deep artificial intelligence structure for online soft measurement that incorporates a nonlinear output error (NOE) framework and a nonlinear autoregressive with an exogenous input (NARX) predictor. This structure provides key insights into the primary characteristics of simulated moving bed chromatography equipment [22]. Santos et al. have developed a theoretical model to describe the adsorption and desorption kinetics of succinic acid and water mixtures in commercial resin-fixed beds. This model has been utilized for the design and optimization of simulated moving bed control processes [23]. Suzuki et al. employed parameter estimation methods to estimate the system, validated the predictive capability of the resulting model using test datasets, and evaluated the confidence intervals of the parameters. These tests confirmed that the model shows improvement compared to the initially obtained model [24]. Other relevant studies can be referenced in [25,26,27].

Overall, past research can be categorized into two main aspects: the study of computational models for SMB systems and the research on controllers. Most computational models focus on the finite element method (FEM), as FEM theory is well established, and many controller designs are based on the FEM. However, this raises two significant issues. Although the FEM provides high accuracy, it suffers from low computational efficiency, leading to prolonged latency in controllers and affecting their performance. Moreover, FEM-based controllers cannot achieve real-time online control with existing hardware infrastructure. To address these issues, this study initially explores the application of the Crank–Nicolson discretization method in SMB systems, demonstrating that the finite difference method (FDM) is also feasible. Subsequent chapters of the study mathematically validate the convergence of iterative calculations and the order of errors. Experimental simulations are then conducted, adjusting parameters of the SMB system, including feed concentration, adsorbent parameters, and switching times, to verify the feasibility of the finite difference method. Finally, a PID controller is used to control the SMB purity based on the finite difference method, further confirming that controllers based on finite difference methods are also stable and reliable.

2. Discrete Model for SMB

2.1. Using Crank–Nicolson Method for Discrete Partial Differential Equations

The SMB model originates from the TMB model. The meanings of the parameters are shown in Table 1 below.

The relevant literature [10] provides a brief introduction to the transformation between the TMB model and the SMB model.

(1) $\frac{\partial C_{i j}}{\partial t} = D_{i} \frac{\partial^{2} C_{i j}}{\partial x^{2}} - v_{j}^{T M B} \frac{\partial C_{i j}}{\partial x} - \frac{1 - ε}{ε} k_{i} (q_{i j}^{*} - q_{i j})$

(2) $\frac{\partial q_{i j}}{\partial t} = \frac{\partial}{\partial x} u_{s} q_{i j} + k_{i} (q_{i j}^{*} - q_{i j})$

$v_{j}^{T M B}$ is the velocity of the TMB system. Therefore, based on the principle of kinematic synthesis, the formula for the conversion between velocities is expressed as Equation (3).

(3) $v_{j}^{S M B} = v_{j}^{T M B} + u_{s}$

Since the solid phase is achieved through column switching controlled by switching time, the equivalent flow rate of the solid phase can be expressed by Equation (4).

(4) $u_{s} = \frac{L}{T_{θ}}$

$L$ is the length of the column; the mathematical model of the SMB system can be obtained as:

(5) $\frac{\partial C_{i j}}{\partial t} = D_{i} \frac{\partial^{2} C_{i j}}{\partial x^{2}} - v_{j}^{S M B} \frac{\partial C_{i j}}{\partial x} - \frac{1 - ε}{ε} k_{i} (q_{i j}^{*} - q_{i j})$

(6) $\frac{\partial q_{i j}}{\partial t} = k_{i} (q_{i j}^{*} - q_{i j})$

$v_{j}^{S M B}$ is denoted as $v_{j}^{*}$ , and we can obtain Equation (7) by substituting Equation (6) into Equation (5).

(7) $\frac{\partial C_{i j}}{\partial t} = D_{i} \frac{\partial^{2} C_{i j}}{\partial x^{2}} - v_{j}^{*} \frac{\partial C_{i j}}{\partial x} - \frac{1 - ε}{ε} \frac{\partial q_{i j}}{\partial t}$

In this study, the focus is mainly on the case of linear isotherms, as follows:

(8) $q_{i j} = H_{i} C_{i j}$

Set

(9) $E = \frac{1 - ε}{ε}$

Equation (7) can be rewritten by substituting Equations (8) and (9), shown as Equation (10).

(10) $(1 + E H_{i}) \frac{\partial C_{i j}}{\partial t} = D_{i} \frac{\partial^{2} C_{i j}}{\partial x^{2}} - v_{j}^{*} \frac{\partial C_{i j}}{\partial x}$

The Crank–Nicolson method employs central differencing for the spatial domain and utilizes the trapezoidal rule for the temporal domain, thereby achieving second-order convergence in time. For instance, for a one-dimensional partial differential equation,

(11) $\frac{\partial u}{\partial t} = F (u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}})$

Denote $Δ x$ and $Δ t$ as the spatial and temporal step sizes, respectively, and denote $u (i Δ x, j Δ t) = u_{i}^{j}$ , then the Crank–Nicolson discretized formula is as follows [28]:

(12) $\frac{u_{i}^{j + 1} - u_{i}^{j}}{Δ t} = \frac{1}{2} (F_{i}^{j + 1} (u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}}) + F_{i}^{j} (u, x, t, \frac{\partial u}{\partial x}, \frac{\partial^{2} u}{\partial x^{2}}))$

Since the Crank–Nicolson method uses central difference formulas in the spatial domain, its formulas are expressed as follows:

(13) $\frac{\partial^{2} u}{\partial x^{2}} |_{i} = \frac{u_{i + 1} - 2 u_{i} + u_{i - 1}}{Δ x^{2}}$

(14) $\frac{\partial u}{\partial x} |_{i} = \frac{u_{i + 1} - u_{i - 1}}{2 Δ x}$

Comparing Equations (10) and (11), it can be observed that

(15) $(1 + E H_{i}) \frac{\partial C_{i j}}{\partial t} = F (C_{i j}, x, t, \frac{\partial C_{i j}}{\partial x}, \frac{\partial^{2} C_{i j}}{\partial x^{2}})$

(16) $F (C_{i j}, x, t, \frac{\partial C_{i j}}{\partial x}, \frac{\partial^{2} C_{i j}}{\partial x^{2}}) = D_{i} \frac{\partial^{2} C_{i j}}{\partial x^{2}} - v_{j}^{*} \frac{\partial C_{i j}}{\partial x}$

The form of Equation (16) now includes first-order partial derivatives in the spatial domain. Since indices $i, j$ have already been used to represent material types and column numbers, we chose $l, k$ as indices for the spatial and temporal domains, respectively. Therefore, the left-hand side of Equation (15) can be written as

$(1 + E H_{i}) \frac{C_{i j, l}^{k + 1} - C_{i j, l}^{k}}{Δ t}$

Hence, we can have the discretization formula for the right-hand side of Equation (16) as

(17) $\begin{array}{l} \frac{1}{2} (F_{l}^{k + 1} (C_{i j}, x, t, \frac{\partial C_{i j}}{\partial x}, \frac{\partial^{2} C_{i j}}{\partial x^{2}}) + F_{l}^{k} (C_{i j}, x, t, \frac{\partial C_{i j}}{\partial x}, \frac{\partial^{2} C_{i j}}{\partial x^{2}})) \\ = \frac{1}{2} D_{i} (\frac{\partial^{2} C_{i j}}{\partial x^{2}} |_{l}^{k + 1} + \frac{\partial^{2} C_{i j}}{\partial x^{2}} |_{l}^{k}) - \frac{1}{2} v_{j}^{*} (\frac{\partial C_{i j}}{\partial x} |_{l}^{k + 1} + \frac{\partial C_{i j}}{\partial x} |_{l}^{k}) \end{array}$

According to Equations (13) and (14), we can obtain

(18) $\frac{1}{2} D_{i} (\frac{\partial^{2} C_{i j}}{\partial x^{2}} |_{l}^{k + 1} + \frac{\partial^{2} C_{i j}}{\partial x^{2}} |_{l}^{k}) = D_{i} \frac{((C_{i j, l - 1}^{k + 1} - 2 C_{i j, l}^{k + 1} + C_{i j, l + 1}^{k + 1}) + (C_{i j, l - 1}^{k} - 2 C_{i j, l}^{k} + C_{i j, l + 1}^{k}))}{2 Δ x^{2}}$

(19) $\frac{1}{2} v_{j}^{*} (\frac{\partial C_{i j}}{\partial x} |_{l}^{k + 1} + \frac{\partial C_{i j}}{\partial x} |_{l}^{k}) = v_{j}^{*} \frac{((C_{i j, l + 1}^{k + 1} - C_{i j, l - 1}^{k + 1}) + (C_{i j, l + 1}^{k} - C_{i j, l - 1}^{k}))}{4 Δ x}$

By substituting Equations (18) and (19) into Equation (12), we obtain

(20) $\begin{array}{l} (1 + E H_{i}) \frac{C_{i j, l}^{k + 1} - C_{i j, l}^{k}}{Δ t} = D_{i} \frac{((C_{i j, l - 1}^{k + 1} - 2 C_{i j, l}^{k + 1} + C_{i j, l + 1}^{k + 1}) + (C_{i j, l - 1}^{k} - 2 C_{i j, l}^{k} + C_{i j, l + 1}^{k}))}{2 Δ x^{2}} - \\ v_{j}^{*} \frac{((C_{i j, l + 1}^{k + 1} - C_{i j, l - 1}^{k + 1}) + (C_{i j, l + 1}^{k} - C_{i j, l - 1}^{k}))}{4 Δ x} \end{array}$

According to Equation (20), move the terms at time point K + 1 to the left side of the equation and move the terms at time point K to the right side, resulting in the iterative computation Equation (21) as follows:

(21) $\begin{array}{l} (1 + E H_{i} + \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, l}^{k + 1} - (\frac{v_{j}^{*} Δ t}{4 Δ x} + \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, l - 1}^{k + 1} + (\frac{v_{j}^{*} Δ t}{4 Δ x} - \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, l + 1}^{k + 1} \\ = (\frac{v_{j}^{*} Δ t}{4 Δ x} + \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, l - 1}^{k} + (1 + E H_{i} - \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, l}^{k} + (\frac{D_{i} Δ t}{2 Δ x^{2}} - \frac{v_{j}^{*} Δ t}{4 Δ x}) C_{i j, l + 1}^{k} \end{array}$

2.2. Boundary Numerate Condition

With the boundary conditions:

Initial condition:

(22) $C_{i j} (x, 0) = C_{0 i j}$

Space conditions:

(23) $\frac{\partial C_{i j} (x, t)}{\partial x} |_{x = L_{e n d}} = 0$

(24) $D_{i} \frac{\partial C_{i j} (x, t)}{\partial x} |_{x = L_{0}} = v_{j}^{*} [C_{i j} (L_{0}, t) - {\bar{C}}_{i j}^{\sec t} (t)]$

$C_{0 i j}$ represents the initial concentration distribution at time = 0. The space boundaries are represented by $L_{e n d}$ and $L_{0}$ , which are the conditions of the end and initial positions of the pipe column. The ${\bar{C}}_{i j}^{\sec t} (t)$ term is related to the region in which it is located. Since Equation (24) describes the concentration changes at the pipe head, the head of each region depends on the flow rate at the inlet of that region, the material concentration, the concentration at the end of the previous column, and the flow rate of the region in contact with the inlet. Therefore, this term can be divided into three cases, as shown in the following formulas:

(25) ${\bar{C}}_{i j}^{I} (t) = \frac{Q_{I V} C_{i j - 1} (l_{n - 1}, t)}{Q_{I}}, Z o n e I, 1^{s t} c o l u m n$

(26) ${\bar{C}}_{i j}^{III} (t) = \frac{Q_{I I} C_{i j - 1} (l_{n - 1}, t) + Q_{f} C_{f i}}{Q_{I I I}}, Z o n e I I I, 1^{s t} c o l u m n$

(27) ${\bar{C}}_{i j}^{\sec t} (t) = C_{i j - 1} (l_{n - 1}, t), o t h e r$

$C_{f i}$ is the concentration of the material in the imported column, and Q is the flow rate of each section.

Assuming that each column is divided into n + 1 components, Equation (23) can be obtained from Equation (12):

(28) $0 = \frac{\partial C_{i j} (x, t)}{\partial x} |_{x = L_{e n d}} = \frac{1}{2} (\frac{\partial C_{i j}}{\partial x} |_{n}^{k + 1} + \frac{\partial C_{i j}}{\partial x} |_{n}^{k}) = \frac{((C_{i j, n + 1}^{k + 1} - C_{i j, n - 1}^{k + 1}) + (C_{i j, n + 1}^{k} - C_{i j, n - 1}^{k}))}{4 Δ x}$

Through the separation of space points n + 1 and space point n − 1, we obtain

(29) $C_{i j, n + 1}^{k} + C_{i j, n + 1}^{k + 1} = C_{i j, n - 1}^{k} + C_{i j, n - 1}^{k + 1}$

By setting the spatial point $j = n$ in Equation (21), we obtain

(30) $\begin{array}{l} (1 + E H_{i} + \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, n}^{k + 1} - (\frac{v_{j}^{*} Δ t}{4 Δ x} + \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, n - 1}^{k + 1} + (\frac{v_{j}^{*} Δ t}{4 Δ x} - \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, n + 1}^{k + 1} \\ = (\frac{v_{j}^{*} Δ t}{4 Δ x} + \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, n - 1}^{k} + (1 + E H_{i} - \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, n}^{k} + (\frac{D_{i} Δ t}{2 Δ x^{2}} - \frac{v_{j}^{*} Δ t}{4 Δ x}) C_{i j, n + 1}^{k} \end{array}$

Substituting Equation (29) into Equation (30), the following equation can be obtained:

(31) $\begin{array}{l} (1 + E H_{i} + \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, n}^{k + 1} - (\frac{v_{j}^{*} Δ t}{4 Δ x} + \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, n - 1}^{k + 1} + (\frac{v_{j}^{*} Δ t}{4 Δ x} - \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, n - 1}^{k + 1} \\ = (\frac{v_{j}^{*} Δ t}{4 Δ x} + \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, n - 1}^{k} + (1 + E H_{i} - \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, n}^{k} + (\frac{D_{i} Δ t}{2 Δ x^{2}} - \frac{v_{j}^{*} Δ t}{4 Δ x}) C_{i j, n - 1}^{k} \end{array}$

Merging similar items, we obtain

(32) $- \frac{D_{i} Δ t}{Δ x^{2}} C_{i j, n - 1}^{k + 1} + (1 + E H_{i} + \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, n}^{k + 1} = \frac{D_{i} Δ t}{Δ x^{2}} C_{i j, n - 1}^{k} + (1 + E H_{i} - \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, n}^{k}$

Equation (33) can be obtained from Equation (24).

(33) $\begin{array}{l} \frac{v_{j}^{*}}{D_{i}} [C_{i j, 1}^{k} - {\bar{C}}_{i j}^{\sec t} (k)] = \frac{\partial C_{i j} (x, t)}{\partial x} |_{x = L_{0}} = \frac{1}{2} (\frac{\partial C_{i j}}{\partial x} |_{1}^{k + 1} + \frac{\partial C_{i j}}{\partial x} |_{1}^{k}) = \\ \frac{((C_{i j, 2}^{k + 1} - C_{i j, 0}^{k + 1}) + (C_{i j, 2}^{k} - C_{i j, 0}^{k}))}{4 h} \end{array}$

Through the separation of space point 0 and space points 1 and 2, we obtain

(34) $C_{i j, 2}^{k + 1} + C_{i j, 2}^{k} - \frac{4 h v_{j}^{*}}{D_{i}} (C_{i j, 1}^{k} - {\bar{C}}_{i j}^{\sec t} (k)) = C_{i j, 0}^{k + 1} + C_{i j, 0}^{k}$

By setting the spatial point $j = 1$ in Equation (21), it can be obtained that

(35) $\begin{array}{l} (1 + F H_{i} + \frac{D_{i} s}{h^{2}}) C_{i j, 1}^{k + 1} - (\frac{v_{j}^{*} s}{4 h} + \frac{D_{i} s}{2 h^{2}}) C_{i j, 0}^{k + 1} + (\frac{v_{j}^{*} s}{4 h} - \frac{D_{i} s}{2 h^{2}}) C_{i j, 2}^{k + 1} \\ = (\frac{v_{j}^{*} s}{4 h} + \frac{D_{i} s}{2 h^{2}}) C_{i j, 0}^{k} + (1 + F H_{i} - \frac{D_{i} s}{h^{2}}) C_{i j, 1}^{k} + (\frac{D_{i} s}{2 h^{2}} - \frac{v_{j}^{*} s}{4 h}) C_{i j, 2}^{k} \end{array}$

Substituting Equation (34) into Equation (35), we obtain the following equation:

(36) $\begin{array}{l} (1 + E H_{i} + \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, 1}^{k + 1} + (\frac{v_{j}^{*} Δ t}{4 Δ x} - \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, 2}^{k + 1} - (\frac{v_{j}^{*} Δ t}{4 Δ x} + \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, 2}^{k + 1} = \\ - (\frac{v_{j}^{*} Δ t}{4 Δ x} + \frac{D_{i} Δ t}{2 Δ x^{2}}) \frac{4 Δ x v_{j}^{*}}{D_{i}} (C_{i j, 1}^{k} - {\bar{C}}_{i j}^{\sec t} (k)) + (1 + E H_{i} - \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, 1}^{k} + \\ (\frac{D_{i} Δ t}{2 Δ x^{2}} - \frac{v_{j}^{*} Δ t}{4 Δ x}) C_{i j, 2}^{k} + (\frac{v_{j}^{*} Δ t}{4 Δ x} + \frac{D_{i} Δ t}{2 Δ x^{2}}) C_{i j, 2}^{k + 1} \end{array}$

Again, merging similar items can lead to

(37) $\begin{array}{l} (1 + E H_{i} + \frac{D_{i} Δ t}{Δ x^{2}}) C_{i j, 1}^{k + 1} - \frac{D_{i} Δ t}{Δ x^{2}} C_{i j, 2}^{k + 1} = (\frac{v_{j}^{* 2} Δ t}{D_{i}} + \frac{2 v_{j}^{*} Δ t}{Δ x}) {\bar{C}}_{i j}^{\sec t} (k) + \\ (1 + E H_{i} - \frac{D_{i} Δ t}{Δ x^{2}} - \frac{v_{j}^{* 2} Δ t}{D_{i}} - \frac{2 v_{j}^{*} Δ t}{Δ x}) C_{i j, 1}^{k} + \frac{D_{i} Δ t}{Δ x^{2}} C_{i j, 2}^{k} \end{array}$

To simplify symbols, we set

(38) $m = \frac{v_{j}^{*} Δ t}{4 Δ x}$

(39) $p = \frac{D_{i} Δ t}{2 Δ x^{2}}$

Then

(40) $v_{j}^{*} = \frac{4 m Δ x}{Δ t}$

(41) $D_{i} = \frac{2 Δ x^{2} p}{Δ t}$

Substituting Equations (40) and (41) into Equation (37) yields

(42) $\begin{array}{l} (1 + E H_{i} + 2 p) C_{i j, 1}^{k + 1} - 2 p C_{i j, 2}^{k + 1} = (\frac{16 m^{2} Δ x^{2}}{Δ t^{2}} \frac{Δ t Δ t}{2 Δ x^{2} p} + \frac{8 m Δ x Δ t}{Δ t Δ x}) {\bar{C}}_{i j}^{\sec t} (k) + \\ (1 + E H_{i} - 2 p - \frac{16 m^{2} Δ x^{2}}{Δ t^{2}} \frac{Δ t Δ t}{2 Δ x^{2} p} - \frac{8 m Δ x Δ t}{Δ t Δ x}) C_{i j, 1}^{k} + 2 p C_{i j, 2}^{k} \end{array}$

Simplify Equation (42) as follows:

(43) $\begin{array}{l} (1 + E H_{i} + 2 p) C_{i j, 1}^{k + 1} - 2 p C_{i j, 2}^{k + 1} = \frac{8 m (m + p)}{p} {\bar{C}}_{i j}^{\sec t} (k) + \\ (1 + E H_{i} - 2 p - \frac{8 m (m + p)}{p}) C_{i j, 1}^{k} + 2 p C_{i j, 2}^{k} \end{array}$

Substituting Equations (40) and (41) into Equation (21) yields the middle position of a spatial point equation as follows:

(44) $\begin{array}{l} (1 + E H_{i} + 2 p) C_{i j, l}^{k + 1} - (m + p) C_{i j, l - 1}^{k + 1} + (m - p) C_{i j, l + 1}^{k + 1} \\ = (m + p) C_{i j, l - 1}^{k} + (1 + E H_{i} - 2 p) C_{i j, l}^{k} + (p - m) C_{i j, l + 1}^{k} \end{array}$

Substituting Equations (40) and (41) into Equation (32) yields the end position of a spatial point equation.

(45) $- 2 p C_{i j, n - 1}^{k + 1} + (1 + E H_{i} + 2 p) C_{i j, n}^{k + 1} = 2 p C_{i j, n - 1}^{k} + (1 + E H_{i} - 2 p) C_{i j, n}^{k}$

According to Equations (43)–(45), (44), we denote the matric

(46) $A = {[\begin{matrix} 1 + E H_{i} + 2 p & - 2 p & 0 & \dots & 0 \\ - (m + p) & 1 + E H_{i} + 2 p & (m - p) & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & \dots & - (m + p) & 1 + E H_{i} + 2 p & (m - p) \\ 0 & \dots & 0 & - 2 p & 1 + E H_{i} + 2 p \end{matrix}]}_{n \times n}$

(47) $B = {[\begin{matrix} 1 + E H_{i} - 2 p - \frac{8 m (m + p)}{p} & 2 p & 0 & \dots & 0 \\ (m + p) & 1 + E H_{i} - 2 p & - (m - p) & \dots & 0 \\ ⋮ & ⋮ & ⋮ & ⋮ & ⋮ \\ 0 & \dots & (m + p) & 1 + E H_{i} - 2 p & - (m - p) \\ 0 & \dots & 0 & 2 p & 1 + E H_{i} - 2 p \end{matrix}]}_{n \times n}$

Set

(48) $C_{i j}^{k} = {[C_{i j, 1}^{k}, C_{i j, 2}^{k}, \dots, C_{i j, n}^{k}]}_{1 \times n}^{T}$

(49) $w (k) = {(\begin{matrix} \frac{8 m (m + p)}{p} {\bar{C}}_{i j}^{\sec t} (k) & 0 & \dots & 0 & 0 \end{matrix})}_{^{1 \times n}}^{T}$

From Equations (43)–(45), the final iteration equation can be unified into the following matrix equation.

(50) $A C_{i j}^{k + 1} = B C_{i j}^{k} + w (k)$

3. Stability and Convergence Analysis

3.1. Stability Analysis

When discussing the stability of the Crank–Nicolson method applied to SMB systems, it is crucial to consider the sources of error in the discretization process of partial differential equations (PDEs). Typically, errors fall into two categories: truncation errors due to the discretization of derivative approximations and inherent error amplification of the numerical method itself. To estimate truncation errors, the Taylor error formula can be employed, which provides an estimate of the errors introduced by discretizing derivatives. To explore whether errors are amplified or dampened, it is necessary to closely examine the behavior of the finite difference method. Von Neumann stability analysis is commonly used to assess whether errors are amplified in numerical methods. For the SMB discretization process, this method’s stability requires selecting an appropriate step size or time increment to ensure that the amplification factor measured by the Von Neumann stability analysis does not exceed 1. By carefully considering truncation errors and employing an appropriate step size in Von Neumann stability analysis, we can evaluate the stability of the Crank–Nicolson method when applied to SMB systems.

Set $C_{i j}^{*} (k)$ is the exact solution of Equation (50), and set $y_{i j} (k)$ is the approximate solution obtained by calculation and satisfying $A y_{i j} (k + 1) = B y_{i j} (k) + w (k)$ . The difference between them is $e_{i j} (k) = y_{i j} (k) - C_{i j}^{*} (k)$ , and it is satisfied that

(51) $\begin{array}{l} A e_{i j} (k) = A y_{i j} (k) - A C_{i j}^{*} (k) = B y_{i j} (k - 1) + w (k - 1) - (B C_{i j}^{*} (k - 1) + w (k - 1)) \\ = B e_{i j} (k - 1) \end{array}$

If A is nonsingular, we can obtain

(52) $e_{i j} (k) = A^{- 1} B e_{i j} (k - 1)$

The spectral radius of $A^{- 1} B$ must satisfy $ρ (A^{- 1} B) < 1$ to ensure the error $e_{i j} (k)$ is not amplified.

Theorem 1.

If matrix A is strictly diagonally dominant and it satisfies $| a_{i i} | > \sum_{\begin{array}{l} j = 1 \\ j \neq i \end{array}}^{n} | a_{i j} |$ , (i =1, 2, …, n), then matrix A is nonsingular.

Proof.

If matrix A is a singular matrix, then there is a non-zero vector x that satisfies $A x = 0$ . Set $| x_{1} | = \max {| x_{1} |, | x_{2} | \dots | x_{n} |}$ , $s o | x_{1} | \neq 0$ , so it can obtain

(53) $a_{11} x_{1} + a_{12} x_{2} + \dots a_{1 n} x_{n} = 0, a_{11} = - a_{12} \frac{x_{2}}{x_{1}} - \dots a_{1 n} \frac{x_{n}}{x_{1}}$

(54) $| a_{11} | \leq | a_{12} | | \frac{x_{2}}{x_{1}} | + \dots | a_{1 n} | | \frac{x_{n}}{x_{1}} | \leq \sum_{j = 2}^{n} a_{1 j}$

This contradicts the condition, so matrix A is nonsingular. □

Set $m = \frac{v_{j}^{*} Δ t}{4 Δ x} = O (s), p = \frac{D_{i} Δ t}{2 Δ x^{2}} = O (s)$ represents an infinitesimal of the same order. As long as the time step is sufficiently small, matrices A and B exhibit strict diagonal dominance, making them non-singular. Therefore, we can use the Crank–Nicolson method to establish the stability of the iterative process involved in computing the SMB system. By ensuring that the time step is sufficiently small, the strict diagonal dominance of matrices A and B guarantees their non-singularity, thereby enhancing the stability of the iterative process used in the Crank–Nicolson method for computing the SMB system.

Theorem 2.

If the space step ( $Δ x$ ) and time step ( $Δ t$ ) are both sufficiently small, $A C_{i j} (k + 1) = B C_{i j} (k) + w (k)$ is stable.

Proof.

Hypothesis $λ$ is the eigenvalue of $A^{- 1} B$ , and $v$ is the corresponding eigenvector. Select $| | v | |_{\infty} = 1$ , and the dimension of Matrix A is $n \times n$ . So with $| v_{k} | \leq 1, 1 \leq k \leq n$ , choose component index $v_{l} = 1$ . Indicators are divided into three situations:

(1) $2 \leq l \leq n - 1$ . Accordingly, $A^{- 1} B v = λ v \Rightarrow B v = λ A v$ , so we can get in component $l$ :

(55) $\begin{array}{l} (m + p) v_{l - 1} + (1 + E H - 2 p) v_{l} + (p - m) v_{l + 1} = λ [- (m + p) v_{l - 1} + \\ (1 + E H + 2 p) v_{l} + (m - p) v_{l + 1}] \end{array}$

(56) $If | λ | = \frac{| (m + p) v_{l - 1} + (1 + E H - 2 p) + (p - m) v_{l + 1} |}{| [- (m + p) v_{l - 1} + (1 + E H + 2 p) + (m - p) v_{l + 1}] |} < 1$

As $m = \frac{v_{j}^{*} Δ t}{4 Δ x} = O (s), p = \frac{D_{i} Δ t}{2 Δ x^{2}} = O (s),$ so

$\forall ε > 0, \exists δ, w h e n Δ t < δ (Δ x), w e obtain | m | \leq ε, | p | \leq ε$

For any $Δ t$ , select a space step $Δ x$ that is small enough to have $p > m \Rightarrow Δ x < \frac{D}{v^{*}}$ .

So, we obtain

(57) $(m + p) v_{l - 1} + (p - m) v_{l + 1} < 2 n$

For the split line over the fixed point (1, 1)

(58) $(m + p) v_{l - 1} + (p - m) v_{l + 1} = 2 n$

Rectangle $(v_{l - 1}, v_{l + 1}) \in [- 1, 1] \times [- 1, 1]$ (Figure 2) is under the split line, so it satisfies the inequality (57) except $(1, 1)$ shown in Figure 2. Therefore, $| λ | \leq 1$ can be inferred.

Equality is established only if $v_{l - 1} = 1, v_{l + 1} = 1$ .

Now, it has been proven that if the components of the eigenvector satisfy this condition, then a contradiction will arise with the previous component.

(59) $\begin{array}{l} (m + p) v_{l - 2} + (1 + E H - 2 p) v_{l - 1} + (p - m) v_{l} = [- (m + p) v_{l - 2} + \\ (1 + E H - 2 p) v_{l - 1} + (m - p) v_{l}] \end{array}$

Simplify Equation (59), and we obtain $v_{l - 2} = 1$ .

Any component can be achieved by the same principle $v_{k} = 1, k = 1, \dots n$ .

However, for the first component equation

$(1 + E H - 2 p - \frac{8 m (m + p)}{p}) + 2 p = (1 + E H + 2 p) - 2 p \Rightarrow \frac{8 m (m + p)}{p} = 0$

m = \frac{v_{j}^{*} Δ t}{4 Δ x} > 0, p = \frac{D_{i} Δ t}{2 Δ x^{2}} > 0

, so

\frac{8 m (m + p)}{p} > 0

. The statement contradicts the previous equation

\frac{8 m (m + p)}{p} = 0

Thus, $| λ | < 1$ .

(2) $l = 1$ , which means

(60) $(1 + E H - 2 p - \frac{8 m (m + p)}{p}) + 2 p v_{2} = λ [(1 + E H + 2 p) - 2 p v_{2}] \Rightarrow$

(61) $If | λ | = \frac{(1 + E H - 2 p - \frac{8 m (m + p)}{p}) + 2 p v_{2}}{(1 + E H + 2 p) - 2 p v_{2}} < 1$

(62) $v_{2} < 1 + \frac{2 m (m + p)}{p^{2}}$

As $| | v | |_{\infty} = 1$ , so $| λ | < 1$ .

(3) $l = n$ , which means

(63) $1 + E H - 2 p + 2 p v_{n - 1} = λ [(1 + E H + 2 p) - 2 p v_{n - 1}] \Rightarrow$

(64) $If | λ | = \frac{1 + E H - 2 p + 2 p v_{n - 1}}{[(1 + E H + 2 p) - 2 p v_{n - 1}]} < 1$

(65) $v_{n - 1} < 1$

When $v_{n - 1} = 1 \Rightarrow$ $λ = 1$ .

It has been proven that if the components of the feature vector are in this situation, contradictions will occur.

(66) $\begin{array}{l} (m + p) v_{l - 2} + (1 + E H - 2 p) v_{l - 1} + (p - m) v_{l} = [- (m + p) v_{l - 2} + \\ (1 + E H - 2 p) v_{l - 1} + (m - p) v_{l}] \end{array}$

Like the first case, any component can be obtained by the same principle $v_{k} = 1, k = 1, \dots n$ . So $| λ | < 1$ and be proved. □

3.2. Convergence Analysis

Denote

$\begin{array}{l} \frac{\partial C_{i j}}{\partial x} = C_{x}, \frac{\partial^{2} C_{i j}}{\partial x^{2}} = C_{x x},, \frac{\partial^{3} C_{i j}}{\partial x^{3}} = C_{x x x},, \frac{\partial^{4} C_{i j}}{\partial x^{4}} = C_{x x x x}, \frac{\partial C_{i j}}{\partial t} = C_{t}, \frac{\partial^{2} C_{i j}}{\partial t^{2}} = C_{t t} \\ \frac{\partial^{3} C_{i j}}{\partial t^{3}} = C_{t t t}, \frac{\partial^{2} C_{i j}}{\partial x \partial t} = C_{x t}, \frac{\partial^{3} C_{i j}}{\partial x^{2} \partial t} = C_{x x t}, \frac{\partial^{3} C_{i j}}{\partial x \partial t^{2}} = C_{x t t}, \frac{\partial^{4} C_{i j}}{\partial x^{2} \partial t^{2}} = C_{x x t t} \end{array},$

and

C (l, k)

represents the concentration at position

l

in space at time

k

Theorem 3.

The truncation error of SMB discrete dynamical systems is the sum of the squares of the time step size and the spatial step size: $O (s^{2}) + O (h^{2})$ .

Proof:

Using the backward difference formula,

(67) $C_{t} (l, k) = \frac{C (l, k) - C (l, k - 1)}{s} + \frac{1}{2} s C_{t t} (l, k) - \frac{1}{6} s^{2} C_{t t t} (l, k_{1}), k_{1} \in (k - 1, k) s$

$s$ represents the time step. According to the Taylor series expansion formula,

(68) $C_{x x} (l, k - 1) = C_{x x} (l, k) - s C_{x x t} (l, k) - \frac{1}{2} s^{2} C_{x x t t} (l, k_{2}), k_{2} \in (k - 1, k) s$

(69) $C_{x x} (l, k) = C_{x x} (l, k - 1) + s C_{x x t} (l, k) + \frac{1}{2} s^{2} C_{x x t t} (l, k_{2}), k_{2} \in (k - 1, k) s$

By the central difference formula of the second derivative, $h$ represents the spatial step, and the following equations can be obtained.

(70) $C_{x x} (l, k) = \frac{C (l - 1, k) - 2 C (l, k) + C (l + 1, k)}{h^{2}} + \frac{h^{2}}{12} C_{x x x x} (l_{1}, k)$

(71) $C_{x x} (l, k - 1) = \frac{C (l - 1, k - 1) - 2 C (l, k - 1) + C (l + 1, k - 1)}{h^{2}} + \frac{h^{2}}{12} C_{x x x x} (l_{2}, k - 1)$

$l_{1}, l_{2} \in (l - 1, l) h$

According to the central difference formula for the first derivative,

(72) $C_{x} (l, k) = \frac{C (l + 1, k) - C (l - 1, k)}{2 h} - \frac{h^{2}}{12} [C_{x x x} (l_{3}, k) + C_{x x x} (l_{4}, k)]$

(73) $C_{x} (l, k - 1) = \frac{C (l + 1, k - 1) - C (l - 1, k - 1)}{2 h} - \frac{h^{2}}{12} [C_{x x x} (l_{5}, k - 1) + C_{x x x} (l_{6}, k - 1)]$

$l_{3}, l_{4}, l_{5}, l_{6} \in (l - 1, l) s$

For the first-order partial derivative $\frac{\partial C_{i j}}{\partial x}$ , by taking the Taylor series expansion,

(74) $C_{x} (l, k) = C_{x} (l, k - 1) + s C_{x t} (l, k) + \frac{1}{2} s^{2} C_{x t t} (l, k_{3}), k_{3} \in (k - 1, k) s$

The SMB system equation is as follows:

(75) $(1 + E H) \frac{\partial C}{\partial t} = D \frac{\partial^{2} C}{\partial x^{2}} - v \frac{\partial C}{\partial x}$

Substituting Equations (67)–(74) into Equation (75), Equation (76) can be obtained as follows:

(76) $\begin{array}{l} (1 + E H) [\frac{C (l, k) - C (l, k - 1)}{s} + \frac{1}{2} s C_{t t} (l, k) - \frac{1}{6} s^{2} C_{t t t} (l, k_{1})] = \\ \frac{1}{2} D [\frac{C (l - 1, k) - 2 C (l, k) + C (l + 1, k)}{h^{2}} + \frac{h^{2}}{12} C_{x x x x} (l_{1}, k)] + \\ \frac{1}{2} D [s C_{x x t} (l, k) + \frac{1}{2} s^{2} C_{x x t t} (l, k_{2})] + \\ \frac{1}{2} D [\frac{C (l - 1, k - 1) - 2 C (l, k - 1) + C (l + 1, k - 1)}{h^{2}} + \frac{h^{2}}{12} C_{x x x x} (l_{2}, k - 1)] \\ - \frac{1}{2} v^{*} [\frac{C (l + 1, k) - C (l + 1, k)}{2 h} - \frac{h^{2}}{12} (C_{x x x} (l_{3}, k) + C_{x x x} (l_{4}, k))] - \\ \frac{1}{2} v^{*} [s C_{x t} (l, k) + \frac{1}{2} s^{2} C_{x t t} (l, k_{3}) + \frac{C (l - 1, k - 1) - C (l + 1, k - 1)}{2 h} - \\ \frac{h^{2}}{12} (C_{x x x} (l_{5}, k - 1) + C_{x x x} (l_{6}, k - 1))] \end{array}$

Moving the terms in Equation (76) gives the following Equation (77).

(77) $\begin{array}{l} (1 + E H) \frac{C (l, k) - C (l, k - 1)}{s} - \frac{1}{2} D \frac{C (l - 1, k) - 2 C (l, k) + C (l + 1, k)}{h^{2}} - \\ \frac{1}{2} D \frac{C (l - 1, k - 1) - 2 C (l, k - 1) + C (l + 1, k - 1)}{h^{2}} + \frac{1}{2} v^{*} \frac{C (l + 1, k) - C (l + 1, k)}{2 h} + \\ \frac{1}{2} v^{*} \frac{C (l - 1, k - 1) - C (l + 1, k - 1)}{2 h} = - \frac{1}{2} (1 + E H) s C_{t t} (l, k) + \frac{1}{6} (1 + E H) s^{2} C_{t t t} (l, k_{1}) + \\ \frac{h^{2}}{24} D [C_{x x x x} (l_{1}, k) + C_{x x x x} (l_{2}, k - 1)] + \frac{1}{2} D s C_{x x t} (l, k) + \frac{1}{4} D s^{2} C_{x x t t} (l, k_{2}) + \frac{v^{*} h^{2}}{24} [C_{x x x} (l_{3}, k) \\ + C_{x x x} (l_{4}, k) + C_{x x x} (l_{5}, k - 1) + C_{x x x} (l_{6}, k - 1)] - \frac{1}{2} v^{*} [s C_{x t} (l, k) + \frac{1}{2} s^{2} C_{x t t} (l, k_{3})] \end{array}$

In Equation (77), the left side represents the result of the discretization of the SMB system, while the right side contains the residual error terms involving the time step size $s$ and spatial step size $h$ . By rearranging the right side of Equation (77), the following expression, Equation (78), can be obtained.

(78) $\begin{array}{l} - \frac{1}{2} (1 + E H) s C_{t t} (l, k) + \frac{1}{6} (1 + E H) s^{2} C_{t t t} (l, k_{1}) + \frac{h^{2}}{24} D [C_{x x x x} (l_{1}, k) + C_{x x x x} (l_{2}, k - 1)] + \frac{1}{2} D s C_{x x t} (l, k) + \\ \frac{1}{4} D s^{2} C_{x x t t} (l, k_{2}) + \frac{v h^{2}}{12} [C_{x x x} (l_{3}, k) + C_{x x x} (l_{4}, k - 1)] - \frac{1}{2} v [s C_{x t} (l, k) + \frac{1}{2} s^{2} C_{x t t} (l, k_{3})] = \\ = [- \frac{1}{2} (1 + E H) s C_{t t} (l, k) + \frac{1}{2} D s C_{x x t} (l, k) - \frac{1}{2} v s C_{x t} (l, k)] + \\ [\frac{1}{6} (1 + E H) C_{t t t} (l, k_{1}) + \frac{1}{4} D C_{x x t t} (l, k_{2}) - \frac{1}{4} v C_{x t t} (l, k_{3})] s^{2} + \\ [\frac{1}{24} D C_{x x x x} (l_{1}, k) + \frac{1}{24} D C_{x x x x} (l_{2}, k - 1) + \frac{v}{12} C_{x x x} (l_{3}, k) + \frac{v}{12} C_{x x x} (l_{4}, k - 1)] h^{2} \end{array}$

For errors of Equation (78), the first-order expressions about space and time steps include

$- \frac{1}{2} (1 + E H) s C_{t t} (l, k) + \frac{1}{2} D s C_{x x t} (l, k) - \frac{1}{2} v s C_{x t} (l, k)$

By taking $\frac{\partial}{\partial t}$ for two sides of Equation (75), the following equation can be obtained.

(79) $- \frac{1}{2} (1 + E H) s C_{t t} (l, k) + \frac{1}{2} D s C_{x x t} (l, k) - \frac{1}{2} v s C_{x t} (l, k) = 0$

So, the error tail term can be rewritten as:

$[\frac{1}{6} (1 + E H) C_{t t t} (l, k_{1}) + \frac{1}{4} D C_{x x t t} (l, k_{2}) - \frac{1}{4} v C_{x t t} (l, k_{3})] s^{2} + [\frac{D}{24} C_{x x x x} (l_{1}, k) + \frac{D}{24} C_{x x x x} (l_{2}, k - 1) + \frac{v}{12} C_{x x x} (l_{3}, k) + \frac{v}{12} C_{x x x} (l_{4}, k - 1)] h^{2}$

We conclude the error trail of Crank–Nicolson is $O (h^{2}) + O (s^{2})$ , which means $O (Δ x^{2}) + O (Δ t^{2})$ . □

4. Simulation

In the SMB simulation system, we use a 2-2-2-2 column structure, which means there are two columns in each region. The initial parameters are set as shown in Table 2. The time step is $Δ t = 0.1$ s, and there are 50 space points for each column, so $Δ x = \frac{L}{50} = 0.5$ cm.

The experiments were executed utilizing MATLAB R2016a software on a PC with an Intel Core i7-3770K 3.53 GHz processor and 16 GB of RAM which was sourced of Santa Clara, California, USA. The resultant experimental data comprised approximately 30 MB. Figure 3 depicts the separation process of the SMB, while Figure 4 presents the axial concentration distribution curve.

Based on the data presented in Figure 3 and Figure 4, the observed variations in solute concentrations at the two outlets exhibit a pattern consistent with diffusion-driven material separation processes. The axial concentration profiles demonstrate that the separated materials ultimately achieve a relatively stable distribution, closely aligning with the operational behavior of the simulated moving bed (SMB). This observation underscores the practical stability and feasibility of discretizing the SMB system. In subsequent analysis, we will examine the effects of varying system parameters on the overall stability of this discretization approach.

The parameters of the adsorbent are considered to change slowly over time due to the aging of the material. Figure 5 shows how concentration varies with the adsorption rate. Computational simulations not only confirmed the stability of the iterative model but also produced significant results from the SMB model. When the adsorption rates of the two materials become more similar, the separation efficiency at the outlet decreases. From the perspective of axial concentration changes, the concentration curve of material A shifts towards that of material B, causing the two curves to converge and making it impossible to separate the two materials. According to the axial concentration distribution, the curves only shift parallelly without changing shape. Thus, for linear isotherms, the key to successfully separating the two materials lies in finding an adsorbent that differentiates between the adsorption properties of materials A and B. Based on the speed of the axial curve movement, the adsorption parameters show a relatively rapid variation range.

Next, through simulations, we observe the impact of feed concentration and switching time parameters on the final separation. The feed concentration can be controlled and adjusted before entering the feed inlet. Figure 6 shows the axial variation in feed concentration. It can be seen that the separation effect of the feed concentration not only influences the separation curves at the outlet but also significantly affects the shape of the axial separation curves. The curves exhibit a leftward skew effect; higher feed concentrations result in higher peaks, causing the two separation curves to be closer to each other.

Switching time parameters are also classified as slow-varying parameters. Due to inherent errors in electronic and mechanical clocks, significant discrepancies can accumulate over time. As shown in Figure 7, switching time has a substantial impact on separation results. It is not only a critical factor for successful separation but also determines the concentration values of the separation under unchanged conditions for other parameters. The SMB system is highly sensitive to switching time parameters. As switching time increases, the final concentration of the separation gradually decreases. From the perspective of axial concentration curves, the peaks become less sharp and flatten near the peaks.

From Figure 5, Figure 6 and Figure 7, it is evident that the Crank–Nicolson discretization method demonstrates significant stability. The method maintains the computational model’s stability and reliability, irrespective of variations in adsorbent parameters, feed concentration, or switching time. However, the model exhibits sensitivity to changes in switching time, leading to considerable fluctuations in the curve’s overall shape. This sensitivity underscores that the separation performance of the simulated moving bed system is highly dependent on the precise adjustment of the switching time parameter. Hence, the meticulous tuning of this parameter is critical for optimizing the system’s performance. Future research will concentrate on refining these parameters to further enhance the model’s stability and efficiency.

5. Controller Simulation

The successful implementation of a PID controller is contingent upon the precise tuning of its gains. Optimizing these PID gains is crucial for achieving optimal performance. In practical applications, it is necessary to discretize the continuous controller and sample it appropriately. The discrete form of the PID controller can be expressed as follows:

(80) $u (k) = K_{p} e (k) + K_{d} Δ e (k) + K_{i} \sum_{i = 1}^{k} e (i)$

In the PID controller framework, the selected input variables comprise the error, the first-order error difference, and the integral error. The corresponding formula is articulated as follows:

(81) $e_{1} = d e s i r e d B - C_{E, B}$

(82) $e_{2} = d e s i r e d A - C_{R, A}$

(83) $e_{3} = e_{1} + e_{2}$

(84) $Δ e_{1} = e_{1} (k) - e_{1} (k - 1)$

(85) $Δ e_{2} = e_{2} (k) - e_{2} (k - 1)$

(86) $Δ e_{3} = Δ e_{1} + Δ e_{2}$

(87) $\int e_{1} = \sum_{k} e_{1} (k)$

(88) $\int e_{2} = \sum_{k} e_{2} (k)$

(89) $\int e_{3} = \int e_{1} + \int e_{2}$

To gain a deeper understanding of the evolving trends, a periodic averaging process is applied to the data. The formula used for a single-cycle moving average smoothing is as follows:

(90) ${\overline{C}}_{E, B, t} = \frac{\int_{t - T}^{t} C_{E, B, t} d t}{T}$

(91) ${\overline{C}}_{R, A, t} = \frac{\int_{t - T}^{t} C_{R, A, t} d t}{T}$

${\bar{C}}_{E, B, t}$ is the B material purity of the extract port of average periodic, ${\bar{C}}_{R, A, t}$ is the A material purity of the raffinate port of average periodic, and the control structure is shown in the Figure 8 below.

The flow rate control parameters for zones I, II, and III are configured as specified in Table 3. The resulting performance metrics are illustrated in Figure 9, Figure 10, Figure 11 and Figure 12. These figures provide a comprehensive analysis of the PID controller’s efficacy across various control objectives, as evaluated through three distinct case studies.

In the first set of experiments, the switching time parameter was set to 180 s, with the target control purity for material A set at 96% and for material B set at 94%. The actual control purity achieved was 95.44% for material A and 94% for material B. The results of the control are shown in Figure 9.

In the second set of experiments, the switching time parameter was set to 180 s, with the target control purity for both material A and material B set at 90%. The actual control purities achieved were 89.7% for material A and 89.42% for material B. The control results are shown in Figure 10. It can be observed from Figure 10 that a slight oscillation occurred when controlling the purity of material B.

In the third set of experiments, the switching time parameter was set to 178 s, with the target control purity for material A set at 93% and for material B set at 95%. The actual control purities achieved were 92.88% for material A and 95% for material B. The control performance is displayed in Figure 11.

In the fourth set of experiments, the switching time parameter was set to 180 s, with the target control purity for material A set at 94% and for material B set at 96%. The actual control purities achieved were 93.98% for material A and 96% for material B. The control results are shown in Figure 12.

From Figure 9, Figure 10, Figure 11 and Figure 12, it is evident that applying the controller within the discretized simulation system yields exceptional results. This lays a robust foundation for the subsequent optimization of the controller design. With a stable discretized simulation system, it becomes feasible to compare the various advantages and disadvantages of different controllers. Additionally, real-time online control is now achievable. Future research could leverage the discretized simulation system as a surrogate for the real system, implementing model predictive control.

In summary, the stability of the discretized dynamical system provides a superior platform for control simulation and efficiency improvements. Building on this, the stable discretized dynamical system not only enhances the accuracy of control simulations but also facilitates more robust analysis and development of advanced control strategies. The ability to conduct simulations in a controlled, virtual environment allows for the rigorous testing and fine-tuning of control algorithms without the constraints and risks associated with real-world experiments. Consequently, the insights gained from these simulations can lead to more effective and reliable control solutions when implemented in actual systems. This approach underscores the importance of a well-designed discretized model in advancing both theoretical research and practical applications in control system engineering.

6. Conclusions

Unlike methods based on finite elements, this study employs the Crank–Nicolson method to discretize the SMB system and rigorously establishes the reliability and computational convergence rate of the proposed method from a theoretical perspective. Subsequent simulation experiments were conducted to validate these findings. Initially, the impact of various parameter changes on the discretization process was examined, revealing that the discretized dynamical system of the SMB also exhibits high stability and reliability, with the computational model closely aligning with the actual system. Furthermore, the implementation of a controller within this discretized dynamical system has proven to be entirely feasible. This result provides a solid foundation for future controller optimization and real-time online control applications. The theoretical and practical significance of this study lies in the fact that the stability and reliability of the discretized model highlight its potential for robust control applications and efficient system design. The successful integration of the controller within this framework not only confirms its practicality but also paves the way for advancements in control strategies and performance optimization.

Subsequent research can focus on two aspects. One aspect is the design and development of controllers, specifically real-time online controllers based on finite difference methods and comparing the advantages and disadvantages of various controllers. The other aspect is the optimization of control based on the economic performance of the SMB system, utilizing high-efficiency simulations based on finite difference methods to enhance the feasibility of optimization control.

Author Contributions

Conceptualization, R.-C.H. and C.-F.X.; methodology, C.-F.X.; software, C.-F.X.; validation, R.-C.H. and C.-F.X.; formal analysis, R.-C.H.; resources, R.-C.H.; data curation, C.-F.X. and H.Z.; writing—original draft preparation, C.-F.X.; writing—review and editing, R.-C.H.; visualization, H.Z. and C.-F.X.; supervision, R.-C.H.; project administration, R.-C.H. All authors have read and agreed to the published version of the manuscript.

Data Availability Statement

To request the corresponding research paper data for experimental simulation, please submit your application via the following email address: [email protected].

Conflicts of Interest

The authors declare no conflicts of interest.

Footnotes

Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Figures and Tables

Figure 1. Simulated moving bed operation process.

Figure 2. Space split line of [Forumla omitted. See PDF.].

Figure 3. SMB separation process on extract and raffinate.

Figure 4. Axial concentration distribution curve.

Figure 5. The concentration of A and B versus different adsorption rates ([Forumla omitted. See PDF.]).

Figure 6. The concentration of A and B versus different feed port concentrations ([Forumla omitted. See PDF.]).

Figure 7. The concentration of A and B versus different switching time ([Forumla omitted. See PDF.]).

Figure 8. PID controller for SMB.

Figure 9. Group 1 control results (desired B = 94%, desired A = 96%, switch time = 180 s).

Figure 10. Group 2 control results (desired B = 90%, desired A = 90%, switch time = 180 s).

Figure 11. Group 3 control results (desired B = 95%, desired A = 93%, switch time = 178 s).

Figure 12. Group 4 control results (desired B = 96%, desired A = 94%, switch time = 180 s).

Table 1

Parameters of SMB system.

Parameter	Nomenclature	Parameter	Nomenclature
$x (cm)$	Axial distance	$Q ({cm}^{3} \min^{- 1})$	Volume flow rate
$k ({gL}^{- 1})$	Comprehensive mass transfer constant	$t (s)$	Time
$v^{*} (cm \min^{- 1})$	Effect velocity of body	$D ({cm}^{2} \min^{- 1})$	Effective dispersion coefficient
$u_{s} (cm \min^{- 1})$	Solid flow rate	$ε$	Bulk void fraction
$C ({gL}^{- 1})$	Mobile phase concentration	$i$	Material index: A or B
$q ({gL}^{- 1})$	Solid phase concentration	$j$	Column number: 1, 2, 3, 4, 5, 6, 7, 8
$q^{*} ({gL}^{- 1})$	Solid phase concentration at equilibrium between solid phase and mobile phase

Table 2

The initial parameters of SMB.

Parameter	Value	Parameter	Value
$L (cm)$	25	$C_{f, i} ({gL}^{- 1})$	5
$d (cm)$	0.46	$θ (\min)$	3
$H_{A}$	0.001	$Q_{I} ({cm}^{3} \min^{- 1})$	6.75
$H_{B}$	0.45	$Q_{I I} ({cm}^{3} \min^{- 1})$	6.6
$D_{A} ({cm}^{2} \min^{- 1})$	0.2	$Q_{I I I} ({cm}^{3} \min^{- 1})$	7
$D_{B} ({cm}^{2} \min^{- 1})$	1.265	$Q_{I V} ({cm}^{3} \min^{- 1})$	2
$ε$	0.8	spatial number	50

Table 3

PID control parameter setting.

Region	K_p	K_d	K_i
I	0.8	0.02	0.02
II	0.000097	0.00005	0.00005
III	0.02	0.003	0.003

References

1. Rajendran, A.; Paredes, G.; Mazzotti, M. Simulated moving bed chromatography for the separation of enantiomers. J. Chromatogr. A; 2009; 1216, pp. 709-738. [DOI: https://dx.doi.org/10.1016/j.chroma.2008.10.075]

2. Chin, C.Y.; Wang, N.L. Simulated moving bed equipment designs. Sep. Purif. Rev.; 2007; 33, pp. 77-155. [DOI: https://dx.doi.org/10.1081/SPM-200042081]

3. Kim, J.I.; Wankat, P.C.; Sungyong, M.; Yoon, K. Analysis of “focusing” effect in four zone SMB (simulated moving bed) unit for separation of xylose and glucose from biomass hydrolysate. J. Biosci. Bioeng.; 2009; 108, pp. 65-76. [DOI: https://dx.doi.org/10.1016/j.jbiosc.2009.08.194]

4. Suvarov, P.; Kienle, A.; Nobre, C.; Weireld, G.D.; Wouwer, A.V. Cycle to cycle adaptive control of simulated moving bed chromatographic separation processes. J. Process Control; 2014; 24, pp. 357-367. [DOI: https://dx.doi.org/10.1016/j.jprocont.2013.11.001]

5. Choi, J.H.; Kang, M.S.; Lee, C.G.; Wang, N.H.L.; Mun, S. Design of simulated moving bed for separation of fumaric acid with a little fronting phenomenon. J. Chromatogr. A; 2017; 1491, pp. 75-86. [DOI: https://dx.doi.org/10.1016/j.chroma.2017.02.040]

6. Supelano, R.C.; Barreto, A.G., Jr.; Neto, A.S.A.; Secch, A.R.I. One-step optimization strategy in the simulated moving bed process with asynchronous movement of ports: A VariCol case study. J. Chromatogr. A; 2020; 1634, pp. 1672-1718.

7. Reinaldo, C.S.; Gomes, B.A.; Resende, S.A. Optimal performance comparison of the simulated moving bed process variants based on the modulation of the length of zones and the feed concentration. J. Chromatogr. A; 2021; 1651, 462280.

8. Sá Gomes, P.M.D. Advances in Simulated Moving Bed: New Operating Modes: New Design Methodologies and Product (FlexSMB-LSRE) Development. Ph.D. Thesis; University of Porto (FEUP): Porto, Portugal, 2009.

9. Dunnebier, G.; Weirich, I.; Klatt, K.U. Computationally efficient dynamic modelling and simulation of simulated moving bed chromatographic processes with linear isotherms. Chem. Eng. Sci.; 1998; 53, pp. 2537-2546. [DOI: https://dx.doi.org/10.1016/S0009-2509(98)00076-1]

10. Neto, A.S.A.; Secchi, A.R.; Souza, M.B.; Barreto, A.G. Nonlinear model predictive control applied to the separation of praziquantel in simulated moving bed chromatography. J. Chromatogr. A; 2016; 1470, pp. 42-49. [DOI: https://dx.doi.org/10.1016/j.chroma.2016.09.070]

11. Muhammed, T.; Tokay, B.; Conradie, A. Raising the Research Octane Number using an optimized Simulated Moving Bed technology towards greater sustainability and economic return. Fuel; 2023; 337, 126864. [DOI: https://dx.doi.org/10.1016/j.fuel.2022.126864]

12. Yan, Z.; Wang, J.S.; Wang, S.Y.; Li, S.J.; Wang, D.; Sun, W.Z. Model Predictive Control Method of Simulated Moving Bed Chromatographic Separation Process Based on Subspace System Identification. Math. Probl. Eng.; 2019; 2019, 2391891. [DOI: https://dx.doi.org/10.1155/2019/2391891]

13. Leao, C.P.; Rodrigues, A.E. Transient and steady-state models for simulated moving bed processes: Numerical solutions. Comput. Chem. Eng.; 2004; 28, pp. 1725-1741. [DOI: https://dx.doi.org/10.1016/j.compchemeng.2004.01.007]

14. Majeed, H.; Hillestad, M.; Knuutila, H.; Svendsen, H.F. Predicting aerosol size distribution development in absorption columns. Chem. Eng. Sci.; 2018; 192, pp. 25-33. [DOI: https://dx.doi.org/10.1016/j.ces.2018.07.004]

15. Kim, Y.; Kim, T.; Park, C.; Lee, J.; Cho, H.; Kim, M.; Moon, I. Development of novel flow distribution apparatus for simulated moving bed to improve degree of mixing. Comput. Chem. Eng.; 2022; 156, 107553. [DOI: https://dx.doi.org/10.1016/j.compchemeng.2021.107553]

16. Lee, W.S.; Lee, C.H. Dynamic modeling and machine learning of commercial-scale simulated moving bed chromatography for application to multi-component normal paraffin separation NSTL. Sep. Purif. Technol.; 2022; 288, 120597. [DOI: https://dx.doi.org/10.1016/j.seppur.2022.120597]

17. Li, S.; Wei, D.; Wang, J.S.; Yan, Z.; Wang, S.Y. Predictive control method of simulated moving bed chromatographic separation process based on piecewise affine. Int. J. Appl. Math.; 2020; 50, pp. 1-12.

18. Hoon, O.T.; Woo, K.J.; Hwan, S.S.; Hosoo, K.; Kyungmoo, L.; Min, L.J. Automatic control of simulated moving bed process with deep Q-network. J. Chromatogr. A; 2021; 1647, 462073.

19. Natarajan, S.; Lee, J.H. Repetitive model predictive control applied to a simulated moving bed chromatography system. Comput. Chem. Eng.; 2000; 24, pp. 1127-1133. [DOI: https://dx.doi.org/10.1016/S0098-1354(00)00493-2]

20. Klatt, K.U.; Hanisch, F.; Dunnebier, G. Mode-based control of a simulated moving bed chromatographic process for the separation of frutose and glucose. J. Process Control; 2002; 12, pp. 203-219. [DOI: https://dx.doi.org/10.1016/S0959-1524(01)00005-1]

21. Carlos, V.; Alain, V.W. Combination of multi-model predictive control and the wave theory for the control of simulated moving bed plants. J. Chem. Eng. Sci.; 2017; 66, pp. 632-641.

22. Marrocos, H.; Iwakiri, I.G.I.; Martins, M.A.F.; Rodrigues, A.E.; Loureiro, J.M.; Ribeiro, A.M.; Nogueira, I.B.R. A long short-term memory based Quasi-Virtual Analyzer for dynamic real-time soft sensing of a Simulated Moving Bed unit. Appl. Soft Comput.; 2022; 116, 108318. [DOI: https://dx.doi.org/10.1016/j.asoc.2021.108318]

23. Santos, L.U.; Delgado, J.A.; Agueda, V.I. Recovery of a Succinic, Formic, and Acetic Acid Mixture from a Model Fermentation Broth by Simulated Moving Bed Adsorption with Methanol as a Desorbent. Ind. Eng. Chem. Res.; 2022; 61, pp. 672-683. [DOI: https://dx.doi.org/10.1021/acs.iecr.1c03388]

24. Suzuki, K.; Harada, H.; Sato, K.; Okada, K.; Tsuruta, M.; Yajima, T.; Kawajiri, Y. Utilization of operation data for parameter estimation of simulated moving bed chromatography. J. Adv. Manuf. Process.; 2022; 4, 10103. [DOI: https://dx.doi.org/10.1002/amp2.10103]

25. Yang, Y.; Chen, X.; Zhang, N. Optimizing control of adsorption separation processes based on the improved moving asymptotes algorithm. Adsorpt. Sci. Technol.; 2018; 36, pp. 1716-1733. [DOI: https://dx.doi.org/10.1177/0263617418804001]

26. Nogueira, I.B.; Martins, M.A.; Rodrigues, A.E.; Loureiro, J.M.; Ribeiro, A.M. Novel Switch Stabilizing Model Predictive Control Strategy Applied in the Control of a Simulated Moving Bed for the Separation of Bi-Naphthol Enantiomers. Ind. Eng. Chem. Res.; 2020; 59, pp. 1979-1988. [DOI: https://dx.doi.org/10.1021/acs.iecr.9b05238]

27. Xie, C.-F.; Hong, Z.; Rey-Chue, H. Discrete Dynamic System Modeling for Simulated Moving Bed Processes. Mathematics; 2024; 12, 1520. [DOI: https://dx.doi.org/10.3390/math12101520]

28. Crank, J.; Nicolson, P. A Practical Method for Numerical Evaluation of Solutions of Partial Differential Equations for the Heat Conduction Type; Cambridge University Press: Cambridge, UK, 1947; Volume 43, pp. 50-67.

© 2024 by the authors. Licensee MDPI, Basel, Switzerland. This article is an open access article distributed under the terms and conditions of the Creative Commons Attribution (CC BY) license (https://creativecommons.org/licenses/by/4.0/). Notwithstanding the ProQuest Terms and Conditions, you may use this content in accordance with the terms of the License.

Stability and Convergence Analysis of the Discrete Dynamical System for Simulating a Moving Bed

Content area