A Finite Difference Method Using High-Order Schemes to Simulate an Equilibrium-Dispersive Model of Non-Linear Chromatography

Abstract

High-performance liquid chromatography (HPLC) is a dynamic separation process with a lot of parameters having different roles. The equilibrium-dispersive model is relevant for simulating HPLC because it is relatively simple and suitable for high-efficiency processes. The partial differential equation was simulated in many different methods such as semi-analytical methods, finite element methods, and finite difference methods. Many studies using finite difference methods have used the first-order and second-order schemes, but higher-order schemes have not been reported yet. This work is about solving the equation of the equilibrium-dispersive model, using a finite difference method with high-order schemes. The fourth-order central difference scheme was used for estimating diffusion and the fifth-order upwind schemes were used for simulating advection. The model was evaluated by assessing the area recovery of the peak, testing the non-retained substance behavior, and comparing the calculation results with the experimental data. The solutions of the equation will indicate the effects of the operation parameters on the system suitability ones and can be used to predict the behavior of an HPLC system and calculate the system suitability parameters of a novel method set.

1. Introduction

HPLC is a popular analytical method in the pharmaceutical, cosmetic, food, and beverage industries. To develop a novel HPLC method that satisfies a given set of criteria, the researcher must do a lot of experiments and adjust many operation parameters. Modeling and simulation are powerful tools to reduce the number of experiments by predicting the experimental outcome. There are many ways to simulate HPLC; two main methodologies are data-driven methods and theory-based ones. The data-driven methods need a massive amount of experimental data. In the studies [1,2,3,4,5], regression tools such as artificial neural networks (ANN), adaptive neuro-fuzzy inference systems (ANFIS), quantitative structure-retention relationships (QSRR), and multiple linear regression (MLR) were employed to interpret the data. In the age of continuously progressing computer science, these data-based methods will have a huge advantage. They can even simulate the gradient elution mode [2,6] which is extremely complicated when using the other method. However, not everyone can perform a lot of experiments or have access to a comprehensive database. Furthermore, most of the studies can only predict the retention behavior of the analytes, and some of them can calculate the peak width at best. The other approach is using theory and laws to set up equations that describe the operation of the equipment. These models, if well-established, will be very powerful at predicting results and only need a few experimental data to calibrate some parameters. The most relevant among these models are the rate model, plate model, and equilibrium-dispersive model [7]. Plate models are based on the idea that the column consists of many theoretical plates stacked upon each other [8,9,10,11]. This theory can be used to establish a set of n equations or an iterative algorithm to calculate the concentration of the analyte in the mobile phase and stationary phase. However, most of the plate models are simple and cannot simulate competitive adsorption when injecting a mixture of analytes [12]. The rate models are the comprehensive models for simulating chromatography, which take into account the effect of the mass transfer resistance. Each type of rate model has a different degree of complexity [12]. The rate models are often expressed by two equations: one demonstrates the concentration of the analyte in the mobile phase and the other describes the mass transfer between the mobile phase and the surface of the stationary phase. On the other hand, equilibrium-dispersive models are the simplified version of rate models, which neglect the effect of the mass transfer resistance [13,14]. The equilibrium-dispersive models are relevant for simulating the high-efficiency column with minimal mass transfer resistance. In this work, equilibrium is chosen for its simplicity and the high-efficiency nature of the HPLC process. The numerical methods for simulating partial differential equations can be divided into two major types: the finite difference method, and the finite element method. The finite element methods derive numerical solutions for the partial differential equation by calculating the integration [12]. These methods work the best in the case of sudden changes, which are common in mechanical problems. However, the finite element methods do not show sufficient accuracy in solving the problems that feature fluid flow. On the other hand, the finite difference methods (FDM) calculate the derivative, so they are more suitable for cases with gradual changes. However, there are not many studies reporting the application of FDM in chromatography modeling. Only few studies used FDM to simulate the partial differential equation with various degrees of success [15,16,17,18]. In all these works, only the first-order scheme is applied. Because of the high dispersion of the first-order scheme, the diffusion term has to be removed. The increments in time and space have to be determined from the diffusion coefficient so that the diffusion is approximated by the dispersion of the numerical scheme. This approach has several drawbacks. First, the physical meaning is not truly reflected by the numerical method. Second, the shape of the peaks sometimes cannot be precisely reproduced. For example, the study from [16] had to adjust the time and space steps by fitting with the result from other methods to obtain a better peak shape. Third, this method is challenging to extend to a multicomponent system because of the variation of the diffusion coefficient. Moreover, the accuracy of these low-order schemes is quite low. This shortcoming was offset by using extremely small temporal and spatial increments, which dramatically prolonged the calculation time. A potential approach to address the issues is to use a higher-order scheme, especially with the constant increase of computer power. Although using a high-order scheme to solve PDE has existed for many years, the application in chromatographic processes is not yet reported. The applicability of the FDM with a high-order scheme to reproduce the peak is still a question. Therefore, the investigation of a high-order-scheme FDM to model the chromatographic processes is worth conducting. In this work, fourth-order central difference schemes and fifth-order upwind schemes are employed to simulate the partial differential equation describing non-linear chromatography with high accuracy, exactly reproducing its peak shape and fast calculation time.

2. Modeling and Simulation

2.1. Mathematical Model

This work uses a rate model with two equations; the first one demonstrates the mass balance in the mobile phase [12]:

$$\frac{\partial C}{\partial t}+\theta \frac{\partial n}{\partial t}=-u\frac{\partial C}{\partial x}+D\frac{{\partial }^{2}C}{\partial x^2}$$

(1)

Here, ν = (1 − ϵ)/ϵ is the phase ratio; C and n are the molar liquid phase and adsorbed phase concentration of the compound, respectively; μ is the dimensionless interstitial velocity; and τ and x are the dimensionless time and space variables.

$$\theta =\frac{6\left(1-\epsilon \right)}{d_p\epsilon }$$

(2)

According to equilibrium theory [15,19,20], the concentration of the analyte in the surface of the stationary phase achieves equilibrium with the concentration in the mobile phase immediately. Thus, the concentration in the stationary phase surface only depends on the concentration in the mobile phase. According to the chain rule of the derivative of the composite function, the derivative of the concentration on the stationary phase surface can be calculated by the derivative of the concentration in the mobile phase, by the following equation:

$$\frac{\partial n}{\partial t}=\frac{\partial n}{\partial C}\times \frac{\partial C}{\partial t}$$

(3)

Substituting Equation (3) to Equation (1), a partial differential equation can be obtained, that can be solved by finite difference methods:

$$\left(1+\theta \frac{\partial n}{\partial C}\right)\frac{\partial C}{\partial t}=-u\frac{\partial C}{\partial x}+D\frac{{\partial }^{2}C}{\partial x^2}$$

(4)

Assigning $r=\theta \frac{\partial n}{\partial C}$, we obtain

$$\frac{\partial C}{\partial t}=-\frac{u}{1+r}\frac{\partial C}{\partial x}+\frac{D}{1+r}\frac{{\partial }^{2}C}{\partial x^2}$$

(5)

The original Langmuir isotherm is as follows:

$$n=\frac{n_0{K}_{L}C}{1+K_{L}C}$$

(6)

One way to increase the flexibility of the model is using the Bi-Langmuir isotherm [20]:

$$n=\frac{n_{01}{K}_{L1}C}{1+K_{L1}C}+\frac{n_{02}{K}_{L2}C}{1+K_{L2}C}$$

(7)

The Langmuir isotherm has a solid theoretical foundation. However, to better fit the experimental data, sometime, a modified Langmuir isotherm must be used. The extended form of the Langmuir isotherm is presented in the following equation [20]:

$$n=\lambda C+\frac{n_0{K}_{L}C}{1+K_{L}C}$$

(8)

For the multi-compound separation, when considering the interaction between different analyte molecules, the competitive Langmuir adsorption isotherm [15,20] is used in the form of the following equation:

$$n_i=\frac{n_{0i}{K}_{Li}{C}_i}{1+{\displaystyle {\sum }_{j=1}^m{K}_{Lj}{C}_j}}$$

(9)

The modified isotherm in the case of competitive adsorption can be addressed by the following equation [20]:

$$n_i={\lambda }_i{C}_i+\frac{n_{0i}{K}_{Li}{C}_i}{1+{\displaystyle {\sum }_{j=1}^m{K}_{Lj}{C}_j}}$$

(10)

The porosity of the column is extremely important in the HPLC modeling and simulation. There are several methods to measure the porosity by experimental data [12,20]. One of them is injecting a non-retained compound into the column and using the retention time of this compound to calculate the porosity. In reversed-phase HPLC, thiourea is a popular substance for this purpose. The retention time consists of the time it takes to flow through the guard column, column, and the system of HPLC:

$$R{T}_0=t_{void}+\frac{L_1}{u_1}+\frac{L_2}{u_2}$$

(11)

$$\epsilon =\frac{4F\left(R{T}_0-t_{void}\right)}{\pi \left(I{D}_1{}^2{L}_1+I{D}_2{}^2{L}_2\right)}$$

(12)

2.2. Numerical Scheme

2.2.1. First-Order Derivatives

There are several ways to derive the approximative schemes for first-order derivatives from the Taylor expansion. LeVeque [21] used a method of undetermined coefficients to define the approximations. In this report, an approach was proposed to obtain these coefficients:

First, we start with the general formular of Taylor expansion:

$$C\left(x+\Delta x\right)=C\left(x\right)+\Delta x\frac{\partial C}{\partial x}+\frac{\Delta x^2}{2!}\frac{{\partial }^{2}C}{\partial x^2}+\frac{\Delta x^3}{3!}\frac{{\partial }^{3}C}{\partial x^3}+\frac{\Delta x^4}{4!}\frac{{\partial }^{4}C}{\partial x^4}+\frac{\Delta x^5}{5!}\frac{{\partial }^{5}C}{\partial x^5}+\frac{\Delta x^6}{6!}\frac{{\partial }^{6}C}{\partial x^6}$$

(13)

When multiplying Δx with a factor a, the following can be obtained:

$$C\left(x+a\Delta x\right)=C\left(x\right)+a\Delta x\frac{\partial C}{\partial x}+\frac{a^2\Delta x^2}{2!}\frac{{\partial }^{2}C}{\partial x^2}+\frac{a^3\Delta x^3}{3!}\frac{{\partial }^{3}C}{\partial x^3}+\frac{a^4\Delta x^4}{4!}\frac{{\partial }^{4}C}{\partial x^4}+\frac{a^5\Delta x^5}{5!}\frac{{\partial }^{5}C}{\partial x^5}+\dots$$

(14)

Suppose that we want to gain a formular from the q position of C (x + a_iΔx): we must multiply C (x + a_iΔx) with the respective coefficients b_i, and add these equations together to obtain the following equation:

$$\begin{gathered} {\displaystyle \sum _{i=1}^q{b}_{i}C\left(x+a_i\Delta x\right)}={\displaystyle \sum _{i=1}^q{b}_{i}C\left(x\right)}+{\displaystyle \sum _{i=1}^q{b}_i{a}_i\frac{\partial C}{\partial x}}+{\displaystyle \sum _{i=1}^q{b}_i\frac{a_i{}^2\Delta x^2}{2!}\frac{{\partial }^{2}C}{\partial x^2}}+{\displaystyle \sum _{i=1}^q{b}_i\frac{a_i{}^3\Delta x^3}{3!}\frac{{\partial }^{3}C}{\partial x^3}} \\ \hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}\hspace{1em}+{\displaystyle \sum _{i=1}^q{b}_i\frac{a_i{}^4\Delta x^4}{4!}\frac{{\partial }^{4}C}{\partial x^4}}+{\displaystyle \sum _{i=1}^q{b}_i\frac{a_i{}^5\Delta x^5}{5!}\frac{{\partial }^{5}C}{\partial x^5}+{\displaystyle \sum _{i=1}^q{b}_i\frac{a_i{}^6\Delta x^6}{6!}\frac{{\partial }^{6}C}{\partial x^6}}+\dots } \end{gathered}$$

(15)

Thus,

$$\frac{\partial C}{\partial x}=\frac{{\displaystyle {\sum }_{}^{}{b}_{i}C\left(x+a_i\Delta x\right)}-C\left(x\right){\displaystyle {\sum }_{i=1}^q{b}_i}}{\Delta x{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i}}-\frac{{\displaystyle {\sum }_{m=2}^{\infty }\left(\frac{\Delta x^{m-1}}{m!}\frac{{\partial }^{m}C}{\partial x^m}{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i{}^m}\right)}}{{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i}}$$

(16)

Since we want to approximate a first-order derivative from the value of function C, all the coefficients of the derivatives $\frac{{\partial }^{m}C}{\partial x^m}$ should be zero. Therefore, the approximation is

$$\frac{\partial C}{\partial x}=\frac{{\displaystyle {\sum }_{}^{}{b}_{i}C\left(x+a_i\Delta x\right)}-C\left(x\right){\displaystyle {\sum }_{i=1}^q{b}_i}}{\Delta x{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i}}$$

(17)

or

$$C_x\Delta x=\frac{{\displaystyle {\sum }_{}^{}{b}_{i}C\left(x+a_i\Delta x\right)}-C\left(x\right){\displaystyle {\sum }_{i=1}^q{b}_i}}{{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i}}$$

(18)

However, with q variables, we can eliminate a maximum of q − 1 terms $\frac{{\partial }^{m}C}{\partial x^m}$.The leftover term will be the truncation error:

$$error=\frac{{\displaystyle {\sum }_{m=2+q-1}^{\infty }\left(\frac{\Delta x^{m-1}}{m!}\frac{{\partial }^{m}C}{\partial x^m}{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i{}^m}\right)}}{{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i}}=O\left(\Delta x^q\right)$$

(19)

This means that when q coefficients are used, the order of accuracy is q.

A MATLAB script was written to the calculate the coefficients of the approximations and their respective truncation errors. Using the MATLAB script, we can gain approximation schemes from any order of accuracy, with an arbitrary set of positions.

2.2.2. Second-Order Derivatives

A similar approach could be used to establish the approximation for a second-order derivative:

$$\frac{{\partial }^{2}C}{\partial x^2}=\frac{2{\displaystyle {\sum }_{i=1}^q{b}_{i}C\left(x+a_i\Delta x\right)}-2C\left(x\right){\displaystyle {\sum }_{i=1}^q{b}_i}}{\Delta x^2{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i{}^2}}$$

(20)

or

$$C_{xx}\Delta x^2=\frac{2{\displaystyle {\sum }_{i=1}^q{b}_{i}C\left(x+a_i\Delta x\right)}-2C\left(x\right){\displaystyle {\sum }_{i=1}^q{b}_i}}{{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i{}^2}}$$

(21)

and the error estimation is

$$error=\frac{{\displaystyle {\sum }_{m=3+q-1}^{\infty }\left(\frac{\Delta x^{m-2}}{m!}\frac{{\partial }^{m}C}{\partial x^m}{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i{}^m}\right)}}{{\displaystyle {\sum }_{i=1}^q{b}_i{a}_i{}^2}}$$

(22)

2.2.3. Sample Injection

The sample injection is a high-complexity process. The concentration of the left boundary as a function of time reflects the injection pulse. This model has the following assumptions to simplify the injection model. The concentration signal from the injector to the column and from the column to the detector does not have a significant change. The delayed time of the injection pulse is the system void volume. An injection signal is a rectangular pulse function with a high equal to the injection concentration and length of injection time. The injection time is $t_{inj}=\frac{V_{inj}}{F}$.

2.2.4. Calculation of the Derivative of Concentration by Numerical Method

To apply an adsorption model to the PDE, we must calculate the derivative of concentration of the analyte on the surface of the stationary phase (n) with respect to the concentration of the analyte in the mobile phase (C). Some adsorption models are very complicated, so calculating the derivative can be extremely tricky. This study uses a fourth-order central difference scheme to calculate the derivative with high accuracy:

$$\frac{\partial n}{\partial C}\approx \frac{n_{i-2}-8n_{i-1}+8n_{i+1}-n_{i+2}}{12\Delta C}$$

(23)

Firstly, the approximation scheme must be validated by calculating the derivative of the analyte concentration on the surface of the stationary phase with respect to the analyte concentration in the mobile phase, with regard to the Langmuir model. The validation will provide the optimal value of the concentration increment. Based on the optimal concentration increment, we can choose a value of concentration with which to calculate the derivative. Then, we must calculate a set of concentrations around the concentration that we want to calculate the derivative.

Then, we can calculate the value of n with respect to C, and, finally, calculate the derivative.

2.2.5. Approximation of the Derivative of Concentrations with Respect to Space

The equilibrium-dispersive model is simulated by solving the Equation (4). Defining Courant number $Cr=\frac{u}{\left(1+r\right)}\frac{\Delta t}{\Delta x}=\frac{4F}{\left(1+r\right)\epsilon \pi I{D}^2}\frac{\Delta t}{\Delta x}$ and alpha number $\alpha =\frac{D}{\left(1+r\right)}\frac{\Delta t}{\Delta x^2}$, the following can be obtained:

$$C_t\Delta t=-Cr\hspace{0.17em}{C}_x\Delta x+\alpha \hspace{0.17em}{C}_{xx}\Delta x^2$$

(24)

The First-Order Upwind Scheme

The normal positions will be calculated by the following equation:

$$C_t\Delta t=-Cr\left(-C_{i-1}+C_i\right)+\alpha \left(C_{i-1}-2C_i+C_{i+1}\right)$$

(25)

The right boundary is calculated by the following equation:

$$C_t\Delta t=-Cr\left(-C_{i-1}+C_i\right)+2\alpha \left(C_{i-1}-C_i\right)$$

(26)

The Second-Order Upwind Scheme

The position next to the left boundary is calculated by the following equation:

$$C_t\Delta t=-Cr\left(-C_{i-1}+C_i\right)+\alpha \left(C_{i-1}-2C_i+C_{i+1}\right)$$

(27)

The normal positions will be calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(C_{i-2}-4C_{i-1}+3C_i\right)}{2}+\alpha \left(C_{i-1}-2C_i+C_{i+1}\right)$$

(28)

The right boundary is calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(C_{i-2}-4C_{i-1}+3C_i\right)}{2}+2\alpha \left(C_{i-1}-C_i\right)$$

(29)

The Third-Order Upwind Scheme

The position next to the left boundary is calculated by the following equation:

$$C_t\Delta t=-Cr\left(-C_{i-1}+C_i\right)+\alpha \left(C_{i-1}-2C_i+C_{i+1}\right)$$

(30)

The normal positions will be calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(C_{i-2}-6C_{i-1}+3C_i+2C_{i+1}\right)}{2}+\alpha \left(C_{i-1}-2C_i+C_{i+1}\right)$$

(31)

The right boundary is calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(C_{i-2}-6C_{i-1}+3C_i+2C_{i+1}\right)}{2}+2\alpha \left(C_{i-1}-C_i\right)$$

(32)

The Fourth-Order Upwind Scheme

Position 2 (next to the left boundary) is calculated by the following equation:

$$C_t\Delta t=-Cr\left(-C_{i-1}+C_i\right)+\alpha \frac{\left(8C_{i-1}-15C_i+8C_{i+1}-C_{i+2}\right)}{6}$$

(33)

Position 3 (next to position 2) is calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(C_{i-2}-6C_{i-1}+3C_i+2C_{i+1}\right)}{6}+\alpha \frac{\left(-C_{i-2}+16C_{i-1}-30C_i+16C_{i+1}-C_{i+2}\right)}{12}$$

(34)

The normal positions will be calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(-C_{i-3}+6C_{i-2}-18C_{i-1}+10C_i+3C_{i+1}\right)}{12}+\alpha \frac{\left(-C_{i-2}+16C_{i-1}-30C_i+16C_{i+1}-C_{i+2}\right)}{12}$$

(35)

The position next to the right boundary is calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(-C_{i-3}+6C_{i-2}-18C_{i-1}+10C_i+3C_{i+1}\right)}{12}+\alpha \frac{\left(-C_{i-2}+8C_{i-1}-15C_i+8C_{i+1}\right)}{6}$$

(36)

The right boundary is calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(3C_{i-4}-16C_{i-3}+36C_{i-2}-48C_{i-1}+25C_i\right)}{12}+\alpha \frac{\left(-C_{i-2}+16C_{i-1}-15C_i\right)}{6}$$

(37)

The Fifth-Order Upwind Scheme

Position 2 is calculated by the following equation:

$$C_t\Delta t=-Cr\left(-C_{i-1}+C_i\right)+\alpha \frac{\left(8C_{i-1}-15C_i+8C_{i+1}-C_{i+2}\right)}{6}$$

(38)

Position 3 is calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(C_{i-2}-6C_{i-1}+3C_i+2C_{i+1}\right)}{6}+\alpha \frac{\left(-C_{i-2}+16C_{i-1}-30C_i+16C_{i+1}-C_{i+2}\right)}{12}$$

(39)

The normal position will be calculated by the following equation:

$$C_t\Delta t=-Cr\frac{\left(-2C_{i-3}+15C_{i-2}-60C_{i-1}+20C_i+30C_{i+1}-3C_{i+2}\right)}{60}+\alpha \frac{\left(-C_{i-2}+16C_{i-1}-30C_i+16C_{i+1}-C_{i+2}\right)}{12}$$

(40)

Position next to the right boundary is calculated by following equation:

$$C_t\Delta t=-Cr\frac{\left(-3C_{i-4}-20C_{i-3}+60C_{i-2}-120C_{i-1}+65C_i+12C_{i+1}\right)}{60}+\alpha \frac{\left(-C_{i-2}+8C_{i-1}-15C_i+8C_{i+1}\right)}{6}$$

(41)

Right boundary is calculated by following equation:

$$C_t\Delta t=-Cr\frac{\left(-12C_{i-5}+75C_{i-4}-200C_{i-3}+300C_{i-2}-300C_{i-1}+137C_i\right)}{60}+\alpha \frac{\left(-C_{i-2}+16C_{i-1}-15C_i\right)}{6}$$

(42)

2.2.6. Integration of Concentration with Respect to Time

The approximation of the derivative of concentration with respect to time is calculated by a one-sided scheme. The first row is the initial condition; the concentration of the analyte in the entire column at that moment is zero. The second moment (second time step) is calculated by the following equation:

$$\frac{\partial C}{\partial t}=\frac{-C_{j-1}+C_j}{\Delta t}\to C_j=C_t\Delta t+C_{j-1}$$

(43)

in which C_tΔt is evaluated at the current position by the derivative with respect to the space-variable (x) as discussed in the Section 2.2.

Similarly, the third moment (third time step) is calculated by the following equation:

$$\frac{\partial C}{\partial t}=\frac{C_{j-2}-4C_{j-1}+3C_j}{2\Delta t}\to C_j=\frac{2C_t\Delta t+4C_{j-1}-C_{j-2}}{3}$$

(44)

The fourth moment (fourth time step) is calculated by the following equation:

$$\frac{\partial C}{\partial t}=\frac{-2C_{j-3}+9C_{j-2}-18C_{j-1}+25C_j}{6\Delta t}\to C_j=\frac{6C_t\Delta t+2C_{j-3}-9C_{j-2}+18C_{j-1}}{11}$$

(45)

The fifth moment (fifth time step) and after are calculated by the following equation:

$$\frac{\partial C}{\partial t}=\frac{3C_{j-4}-16C_{j-3}+36C_{j-2}-48C_{j-1}+25C_j}{12\Delta t}\to C_j=\frac{12C_t\Delta t-3C_{j-4}+16C_{j-3}-36C_{j-2}+48C_{j-1}}{25}$$

(46)

2.2.7. Injection Simulation

The sample injection is a high-complexity process. The concentration of the left boundary as a function of time reflects the injection pulse. This model has the following assumptions to simplify the injection model. The concentration signal from the injector to the column and from the column to the detector does not have a significant change. The delayed time of the injection pulse is the system void volume. An injection signal was a Gaussian pulse function presented by the following equation:

$$C\left(t\right)=C_{inj}\exp \left(-\frac{{\left(t-t_{void}\right)}^2}{2\hspace{0.17em}w{d}^2}\right)$$

(47)

The standard deviation of the Gaussian pulse was calculated by the following equation:

$$wd=\frac{V_{inj}}{F\sqrt{2\pi }}$$

(48)

3. Calculation and Experiment

The algorithm for modeling PDEs is written in the MATLAB programming language and run by MATLAB R2021a. The code designed to study the effect of the input parameters on the SST parameters is written in different scripts and run to obtain the results. The parameters which are surveyed are the flow rate, injection volume, injection concentration, stationary phase capacity, adsorption isotherm, column temperature, diffusion coefficient, temporal increment, and spatial increment. The undefined parameters, such as the stationary phase capacity, adsorption isotherm, diffusion coefficient, temporal increment, spatial increment, are calibrated by the experimental data.

Several approximation schemes are derived from the Taylor expansion to simulate the PDE, and their results are put side by side. These schemes are validated by calculating the derivative of some test functions and comparing the numerical results with the analytical ones.

The validation experiments were carried out on the Agilent HPLC system, using the C8 reversed phase column; the column length is 250 mm, inner diameter is 4.6 mm, particle size is 5 µm, detection wavelength is 275 nm, and flow rate is 1 mL/min. The standard solution is caffeine 0.25 mg/mL. Some experiments were performed to calibrate the undefined parameters and validate the model. The HPLC Agilent 1260 Infinity II was used for running the sample. Information about the chemical reagents which were used for the sample preparation is presented in

Table 1. Information about the chemical reagents.

Name	Vendor	Lot/Batch
Methanol (HPLC grade)	J.T.Baker	0000243087
Ortho-phosphoric acid (85%)	Merck	K50782473
Distilled water	HCMC IDQC	NA

.

4. Results and Discussions

4.1. Numerical Case Studies

The distribution of concentration of the analyte in the column at a specific point in time has a belt shape as shown in Figure 1. In the figure, the run time is the time it takes to run an experiment [22], and it is set to be 20% higher than the retention time of the peak in the chromatogram in this study. In the numerical experiment, the retention time of the peak is 9.545 min and the run time is 11.45 min, which is 120% of the retention time to ensure the complete elusion. Figure 1 represents the concentration distribution of the analyte in the column at different moments, namely, 10, 20, 30, 40, 50, 60, 70, 80, and 90% of the run time, corresponding to 1.1450, 2.2900, 3.4350, 4.5800, 5.7250, 6.8700, 8.0150, 9.1600, and 10.3050 min, respectively. They can be either symmetric Gaussian shape peaks when the concentration is in the linear range or slightly tailing when the concentration is in the non-linear range. Over the run time, the analyte shifts toward the column outlet, becomes shorter and broadens, and eventually gets out of the column. The upper chart shows the lower concentration of the injection which is within the linear concentration. In the linear concentration, the distribution of concentration in the column is nearly symmetrical so the peak in the chromatogram is symmetrical as well. By contrast, the lower chart depicts the distribution of concentration in the column above the linear range so the distribution shows significantly more tailing than that on the linear range due to the overloading phenomenon. Therefore, the peak in its chromatogram also shows more tailing.

The comparison of the chromatograms produced by schemes with different orders of accuracy is presented in Figure 2. From the results, a clear improvement can be observed when increasing the order of accuracy. The parameters demonstrating the peak shape and position, namely, the plate count, symmetry factor, and retention time, approach equilibrium values as the order of accuracy goes higher. The area recovery, calculated by dividing the area of the peak with the area of the injection pulse, reflecting the mass conservation properties of the model, is also converted to 100% as the order of accuracy increases. A further inspection of the peak bases depicts an apparent advantage of the higher-order scheme in reducing the oscillation produced by the error of the approximation schemes.

The effect of the absorption isotherm on the chromatogram and suitability parameters are depicted in Figure 3. The retention factor is directly proportional to the absorption isotherm. The retention factor has a linear relationship with the absorption isotherm, and it intercepts the y axis at the retention time of the non-retained component which was used to calculate the porosity of the column. This makes sense because the tracer has an isotherm or/and capacity equal to zero. Running a simulation with either the isotherm or capacity or both equal to zero and comparing the retention time with the dead time is one way to test the model. The increase of the absorption isotherm causes the rise of the asymmetry factor. This is because the absorption isotherm curve will become stiffer and curvier when the absorption isotherm goes up. The higher the isotherm is, the lower the plate count will be. This phenomenon is caused by the combination of two factors, the retention time and asymmetry factor. Firstly, the longer the component retains in the column, the higher the effect of the peak broadening will be. Secondly, the increase of the asymmetry factor is always accompanied by the downfall of the column efficiency and the decrease in plate count.

The effect of the diffusion coefficient on the chromatogram and suitability parameters are shown in Figure 4. From the figure, the diffusion coefficient is responsible for the broadening phenomenon of the peak. When the diffusion coefficient increases, the peak height will drop, and the peak width will broaden until some point where they become relatively stable. This leads to the decrease of the column efficiency and the decline of the plate count. The diffusion coefficient also helps stabilize the peak because the higher the diffusion coefficient, the more symmetric the peak will become. This causes the diffusion coefficient to affect the retention time a little, because the asymmetry factor decrease will slightly raise the retention time.

The effect of the stationary capacity on the retention time and retention factor are shown in Figure 5. Similar to that of the absorption isotherm, the retention factor is directly proportional to the capacity, and the retention factor has a linear relationship with it. By and large, the asymmetry factor gradually decreases when the capacity increases. This is because the longer the component retains in the column, the greater the effect of diffusion will be. Thus, the peak will become more symmetric when the capacity goes up, except for the peak of the unretained component. When the capacity drops, the effect of diffusion decreases, and the peak become more asymmetric. However, the cause of the asymmetry is concentration overload. Concentration overload only occurs when the component interacts with the stationary phase. When the stationary phase capacity is equal to zero, the component no longer interacts with the stationary phase, and, thus, the effect of concentration overload is eliminated. When the stationary phase capacity decreases, the plate also gradually increases. When the capacity becomes zero, the plate counts increase dramatically. This is because the peak suddenly becomes almost symmetrical.

When the injection concentrations are low and vary in a small range as shown in Figure 6, the asymmetry factor, plate count, and retention time do not have any significant changes. However, when the concentration becomes extremely high as presented in Figure 7, it will cause overload, increase the asymmetry factor, decrease the plate count, and shift the peak forward.

The injection volume has a similar effect on the peak shape as the injection concentration. Because HPLC has a limited range of injection volume, compared to preparative chromatography, the volume overload phenomenon does not occur. However, in the linear range, the increase of injection volume makes the peak width a bit larger, and therefore, makes the plate counts decrease a little.

4.2. Model Validation

4.2.1. Assessment of Mass Conservation

To validate the conservative properties of the model, the recovery of the peak area compared to the injection pulse are calculated for each single peak. The area of each peak is calculated by the Composite Trapezoidal Rule. The recovery of peak area is calculated by the following equation:

$$AR\%=\frac{Area\times F}{C_{inj}{V}_{inj}}\times 100\%$$

(49)

The results are shown in

Table 2. Area recovery of the simulation peak.

Area of Injection Pulse	Peak Area	Recovery	Bias
0.7771	0.7765	99.92%	0.08%
1.5541	1.5527	99.91%	0.09%
2.3312	2.3290	99.91%	0.09%
3.1083	3.1054	99.91%	0.09%

. From the results, it can be observed that the recovery of the proposed numerical method is very high, nearly 100%.

4.2.2. Simulation of Non-Retained Substance

In the experiment, the peak of the solvent is 2.73 min, which is used to calculate the column porosity. Then, the porosity is used in the further calculation in this simulation. A chromatogram of the non-retained substance can be simulated by setting all the parameters involved in the interaction between the substance and the stationary phase, such as the stationary phase, absorption isotherm, and Langmuir modified factor, to zeros. The retention time of the simulated tracer is 2.7314 min, which is extremely close to the retention time of the solvent peak. This is evidence that the simulation is operated properly.

4.2.3. Validation by Experiment

Figure 8 presents chromatograms of the experimental and simulation results (from the left to the right). The comparison indicates that the model can simulate a single injection and generate a simulation result which has significant accuracy when comparing the system suitability parameters. From the results in

Table 3. System suitability parameters of the experiment and simulation peak.

Parameters	Experiment	Simulation	Accuracy	Bias
RT, min	9.523	9.543	100.2%	0.2%
As	1.21	1.14	94.5%	5.5%
NEP	14,167	14,283	100.8%	0.8%
NUSP	14,021	14,220	101.4%	1.4%
Wd5, min	0.412	0.3909	94.9%	5.1%
Wd10, min	0.3549	0.3430	96.6%	3.4%
WdUSP, min	0.322	0.3201	99.4%	0.6%

, it can be observed that some parameters such as retention time and tangent peak width have extremely high accuracy, with a bias of not more than 2%. However, this cannot be achieved without compromising other parameters such as asymmetry factor and peak width at 5 and 10%, and EP plate count with a bias of up to 5.5%.

After simulating a single injection with acceptable precision as shown in Table 3, a projection about the effect of some operation parameters can be made and the results are presented in

Table 4. Experiment and simulation result at V_inj = 10 µL.

Parameters	Experiment	Simulation	Accuracy	Bias
RT, min	9.523	9.548	100.3%	0.3%
As	1.21	1.13	93.2%	6.8%
NEP	14,167	14,705	103.8%	3.8%
NUSP	14,021	14,659	104.6%	4.6%
Wd5, min	0.412	0.3851	93.5%	6.5%
Wd10, min	0.3549	0.3379	95.2%	4.8%
WdUSP, min	0.322	0.3155	98.0%	2.0%

,

Table 5. Experiment and simulation result at V_inj = 20 µL.

Parameters	Experiment	Simulation	Accuracy	Bias
RT, min	9.523	9.548	100.3%	0.3%
As	1.21	1.13	93.2%	6.8%
NEP	14,167	14,705	103.8%	3.8%
NUSP	14,021	14,659	104.6%	4.6%
Wd5, min	0.412	0.3851	93.5%	6.5%
Wd10, min	0.3549	0.3379	95.2%	4.8%
WdUSP, min	0.322	0.3155	98.0%	2.0%

,

Table 6. Experiment and simulation result at V_inj = 30 µL.

Parameters	Experiment	Simulation	Accuracy	Bias
RT, min	9.523	9.548	100.3%	0.3%
As	1.21	1.13	93.2%	6.8%
NEP	14,167	14,705	103.8%	3.8%
NUSP	14,021	14,659	104.6%	4.6%
Wd5, min	0.412	0.3851	93.5%	6.5%
Wd10, min	0.3549	0.3379	95.2%	4.8%
WdUSP, min	0.322	0.3155	98.0%	2.0%

and

Table 7. Experiment and simulation result at V_inj = 40 µL.

Parameters	Experiment	Simulation	Accuracy	Bias
RT, min	9.523	9.548	100.3%	0.3%
As	1.21	1.13	93.2%	6.8%
NEP	14,167	14,705	103.8%	3.8%
NUSP	14,021	14,659	104.6%	4.6%
Wd5, min	0.412	0.3851	93.5%	6.5%
Wd10, min	0.3549	0.3379	95.2%	4.8%
WdUSP, min	0.322	0.3155	98.0%	2.0%

. For instance, the effect of injection volume on the system suitability parameters can be forecasted. While the retention times are simulated with errors less than 1%, the accuracy of the peak widths are a bit higher, up to 6.5% of bias at the highest. The increased plate counts when the injection volume decreases are successfully predicted by the model, with a bias of less than 5%. However, the simulation trendlines of the plate count drop more dramatically than that of the experimental peak. The model also predicts the decrease of the asymmetry factor when the injection volume rises, which is quite agreeable with the experiment data. However, the simulation peaks become symmetrical faster than the experimental peaks, and the deviation between the simulation and experiment is quite high, up to 6.5%.

5. Conclusions

The equation representing the analytical process of the HPLC system was solved by a finite difference method with high-order schemes. A procedure to define the coefficients of these schemes was proposed and an algorithm in MATLAB was written to perform said procedure. The algorithm makes it easy to establish schemes at any order with any set of nodes and effortlessly calculate their most significant term of truncation error. The advection and axial diffusion were approximated by the fifth-order upwind schemes and fourth-order central difference scheme, respectively. The derivative with respect to time is approximated by a one-sided fourth-order scheme. The accuracy of these schemes was proven by verification with relevant test functions. Simulation results run by different schemes show the advantage and superiority of the fifth- and fourth-order schemes over their lower-order counterpart.

The effect of the injection volume and injection concentration and other factors on the peak shape were examined by running the model. The simulation results show a linear behavior in a low concentration and low volume of injection. This is represented in highly symmetric peaks, and relatively stable retention time and plate counts. In a high volume and concentration of injection, the system begins to show signs of concentration overload. The symmetry factor of these peaks rapidly increases, and the plate counts swiftly decline. In the multi-component separation, the present of more than one component slightly reduces the capacity of each component. Competitive absorption causes peaks to be shifted forward a little. These behaviors, so far, comply with the well-established and confirmed theories about chromatography and absorption. The model was evaluated by assessing the area recovery of the peak, and the results showed that the simulation performed properly. Experiments for validating the model have been performed. The results show a significant agreement between the experiment and simulation.