2.1 Mathematical model for PSA process
To describe the adsorption process, mass conservation equations, energy conservation equations, momentum equation and adsorption isotherm equation should be written.
The mass conservation for the bulk phase in the adsorption column contains tow parts: for the component and overall mass conservation.
$$-{D}_{L}\frac{{\partial }^{2}{y}_{i}}{\partial {z}^{2}}+\frac{\partial {y}_{i}}{\partial t}+{u}_{z}\frac{\partial {y}_{i}}{\partial z}+\frac{RT}{p}\frac{1-{\epsilon }_{b}}{{\epsilon }_{b}}{\rho }_{p}\left(\frac{\partial {q}_{i}}{\partial t}-{y}_{i}\sum _{j=1}^{N}\frac{\partial {q}_{j}}{\partial t}\right)=0, i=1,···,N$$
1
,
$$-{D}_{L}\frac{{\partial }^{2}p}{\partial {z}^{2}}+\frac{\partial p}{\partial t}+p\frac{\partial u}{\partial z}+{u}_{z}\frac{\partial p}{\partial z}-\frac{p}{T}\left(-{D}_{L}\frac{{\partial }^{2}T}{\partial {z}^{2}}+\frac{\partial T}{\partial t}+{u}_{z}\frac{\partial T}{\partial z}\right)+\frac{1-{\epsilon }_{b}}{{\epsilon }_{b}}{\rho }_{p}\text{R}T\sum _{j=1}^{N}\frac{\partial {q}_{j}}{\partial t}=0$$
2
,
where \({D}_{L}\) is the axial dispersion coefficient; \({u}_{z}\) is the axial physical velocity;\({y}_{i}\) and \({q}_{i}\) represent the molar fraction and the adsorbed phase concentration of species\(i\) respectively; \({\epsilon }_{b}\) is the interparticle void fraction;\({\rho }_{p}\) is the pellet density of the adsorbent;\(\text{R}\) is the universal gas constant, \(T\) and\(P\) correspond to the temperature and pressure in the adsorption bed; \(t\) and \(z\) respectively represent the time since the beginning of the sorption process and the axial position in the bed.
The energy conservation of the PSA system also includes two parts: the energy conservation between the gas and solid phase in the adsorption, the energy conservation at the wall of the adsorption bed.
$$-{K}_{L}\frac{{\partial }^{2}T}{\partial {z}^{2}}+\left({\epsilon }_{t}{\rho }_{g}{C}_{pg}+{\rho }_{b}{C}_{ps}\right)\frac{\partial T}{\partial t}+{\rho }_{g}{C}_{pg}{\epsilon }_{b}{u}_{z}\frac{\partial T}{\partial z}-{\rho }_{b}{\sum }_{i}^{N}{Q}_{i}\frac{\partial {q}_{i}}{\partial t}+\frac{2{ℎ}_{i}}{{R}_{bi}}\left(T-{T}_{w}\right)=0$$
3
$${\rho }_{w}{C}_{pw}{A}_{w}\frac{\partial {T}_{w}}{\partial t}=2\pi {R}_{bi}{ℎ}_{i}\left(T-{T}_{w}\right)-2\pi {R}_{bo}{ℎ}_{o}\left({T}_{w}-{T}_{atm}\right)$$
4
,
$${A}_{w}=\pi \left({R}_{bo}^{2}-{R}_{bi}^{2}\right)$$
5
.
where \({K}_{L}\) is the thermal axial dispersion coefficient;\({C}_{pg}\) is the heat capacity of the gas phase;\({C}_{ps}\) is the specific heat capacity of the adsorbent; \({\epsilon }_{t}\) is the total void fraction of the bed; \({\epsilon }_{b}\) is the interparticle void fraction; \({Q}_{i}\) is the heat of adsorption of species \(i\); \({ℎ}_{i}\) and \({ℎ}_{o}\) are respectively the heat transfer coefficient with inner and outer wall of column; \({T}_{w}\) is the wall temperature and \({T}_{atm}\) is the ambient temperature; \({R}_{bi}\) is the inside radius of adsorption bed; \({R}_{bo}\) is the outside radius of adsorption bed.
Based on the work of Sereno et al. [35], the model neglects the kinetic energy change in the mechanical energy balance. Momentum in the adsorption bed is determined through Ergun’s equation:
$$-\frac{dp}{dz}=a\mu {\upsilon }_{z}+b\rho {\upsilon }_{z}\left|\overrightarrow{\upsilon }\right|$$
6
,
where the coefficients \(\text{a}\) and \(\text{b}\) are determined by the following equations:
$$a=\frac{150}{4{R}_{p}^{2}}\frac{{\left(1-{\epsilon }_{b}\right)}^{2}}{{\epsilon }_{b}^{3}}$$
7
,
$$b=1.75\frac{\left(1-{\epsilon }_{b}\right)}{2{R}_{p}{\epsilon }_{b}^{3}}$$
8
.
In these equations, \(\mu\) is the dynamic viscosity; \({\upsilon }_{z}\) is Darcy’s velocity and \({R}_{p}\) corresponds to the particle radius.
The extended Langmuir-Freundlich model is used in this study to express the adsorption isotherms of multi-component gas:
$${q}_{i}^{\ast }=\frac{{q}_{mi}{B}_{i}{p}_{i}^{{n}_{i}}}{1+{\sum }_{j=1}^{N}{B}_{j}{p}_{j}^{{n}_{j}}}, i=1,\cdots ,N$$
9
.
In this equation, \({q}_{i}^{\ast }\) is the equilibrium absorbed phase concentration;\({q}_{mi} \text{a}\text{n}\text{d} {B}_{i}\) are the extended Langmuir-Freundlich isotherm parameters and\({p}_{i}\) represent the partial pressure of species \(i\). The isotherm parameters are functions of temperature:
$${q}_{m}={k}_{1}+{k}_{2}T, B={k}_{3}{e}^{\frac{{k}_{4}}{T}}$$
10
$$n={k}_{5}+\frac{{k}_{6}}{T}$$
.
2.2 BBD method for PSA process
RSM is a statistical method used to study the relationship between multi factors and response variables, such as Box-Behnken Design (BBD), Central Composite Design (CCD). Within the range of experimental design, the optimal combination of the factors and responses can be found. Under the same factors, the number of BBD experimental design groups is less than that of CCD, making it more economical. BBD is commonly used for predicting nonlinear models [36, 37].
The one factor design was used firstly to determine the range of independent factors for BBD method, each independent optimization parameter is placed at one of three equally spaced values, usually coded as “−1”, “0” and “+1” (the factors). In the BBD method, desirability is an objective function that ranges from zero (outside of the limits of the performance objectives) to one (where the performance objectives are met). A value of one represents the ideal case, and a zero indicates that one or more responses fall outside desirable limits. The overall desirability (D) and the desirability for each response (di) are defined as follows:
$$D={\left({d}_{1}\times {d}_{2}\times \cdots \times {d}_{n}\right)}^{\frac{1}{n}}={\left(\prod _{i=1}^{n}{d}_{i}\right)}^{\frac{1}{n}},$$
11
$${d}_{i}={\left(\frac{{y}_{i}-L}{T-L}\right)}^{w},$$
12
where n is the number of responses in the measure, T and L represent the maximum and minimum possible values for the responses, respectively, w is the weight of each response, and yi is the optimum value of each response as determined by the BBD method.
2.3 BPNN-GA for optimization
ANN is a machine learning method that simulates the structure of the human brain’s nervous system. It is usually used to estimate or approximate functions depend on large number of inputs and generally the functions are unknown. In this work, back propagation neural network (BPNN) is used for prediction and optimization of PSA process, the structure of BPNN is showed in Fig. 1. A tangent sigmoid transfer function (tansig) is used for hidden layer, with a log-sigmoid transfer function (logsig) is used at output layer. For training the designed networks, the Levenberg-Marquardt backpropagation (trainlm) is selected. The calculation formula for forward propagation of BPNN can be expressed as follows:
$${X}_{ij}=\sum _{k=1}^{{N}_{i-1}}{Y}_{\left(i-1\right)k}\times {W}_{\left(i-1\right)kj}$$
13
,
$${Y}_{ij}={f}_{s}\left({net}_{ij}\right)=\frac{1}{1+exp\left[-\left({X}_{ij}-{\theta }_{ij}\right)\right]}$$
14
,
with \({X}_{ij}\) is the total input of the j-th neuron in the i-th layer; \({W}_{ijk}\) is the connection weight from the j-th neuron in the i-th layer to the i + 1-th layer to the k-th neuron; \({Y}_{ij}\) is the output of the j-th neuron in the i-th layer; \({\theta }_{ij}\) is the threshold of the j-th neuron in the i-th layer; \({N}_{i}\) is the number of neurons in the i-th layer.
Although BPNN has strong non-linear mapping ability, it does not guarantee the global optimal solution. Consequently, the genetic algorithm (GA) is introduced to solve the optimization problem in this work, and it is also cost-effective and less time-consuming technique. The procedure of the BPNN optimized by GA is showed in Fig. 2.