L Fangueiro Gomes

F Augier

D Leinekugel-Le-Cocq

I Vinkovic

S Simoëns

Modelling of a Simulated Moving Bed in case of non-ideal hydrodynamics

Keywords: Adsorption, CFD, Simulated Moving Bed, 1D modelling Research Highlights •

. Then, this new model is integrated in a cyclic solver in order to perform Simulated Moving Bed (SMB) studies. The new model provides a detailed hydrodynamic description, which appears to be mandatory especially when mass transfer exchanges are fast, without undergoing the prohibitive simulation times of CFD models.

INTRODUCTION

In adsorption separation processes, the difference of the affinities of an adsorbent for the adsorbates (components present in the feed mixture), generates different retention times, enabling the separation of molecules. However, in the presence of dispersive phenomena, the concentration fronts may be spread, limiting the separation [START_REF] Ruthven | Effect of dead volumes on the performance of an industrial-scale simulated moving-bed Parex unit for p-xylene purification[END_REF]. In some cases, such phenomena can lead to peak tailing: the concentration fronts advance asymmetrically through the separation unit, leaning to the right when plotted in function of time.

The dispersive phenomena mentioned above can essentially be of intra-particle nature and/or hydrodynamic nature. The intra-particle dispersive phenomena comprise the mass transfer resistances between the bulk and adsorbed phases and the non-linearity of the sorption equilibrium. The hydrodynamic dispersive phenomena are related to the adsorption vessel design: packing heterogeneities, injection/withdrawal device and bed geometry.

In order to design and optimize adsorption processes, models are often used to study and quantify the influence of different parameters. One simple approach is to consider that the motion of the eluent can be described by a dispersed plug flow (DPF) model. In this case, the flow is assumed essentially onedimensional and all the hydrodynamic heterogeneities are taken into account by adjusting an axial dispersion coefficient. This approach is sufficiently accurate when the intra-particle dispersive phenomena are dominant relatively to the hydrodynamic ones. However, if intra-particle and hydrodynamic dispersive phenomena are of the same order of magnitude the validity of this model is at stake. In this case a better description of flow may be required. When in the need of a more detailed perception of hydrodynamics of industrial apparatus, it is usual to make use of CFD simulations.

Several studies have already been published regarding the coupling between detailed hydrodynamics and adsorption since the paper of [START_REF] Bart | Two-dimensional non-isothermal model for toluene adsorption in a fixed-bed adsorber[END_REF]. These authors developed a two-dimensional model to study the effect of maldistribution on the sorption of toluene in activated carbon. They found a good agreement between simulation and experimental results. [START_REF] Wu | A computational fluid dynamics study of binary adsorption separation in chromatography[END_REF] and [START_REF] Kwapinski | Thermal and flow effects during adsorption in conventional, diluted and annular packed beds[END_REF] performed CFD simulations in order to evaluate the impact of dispersion on the separation efficiency of a binary separation in a low aspect ratio chromatographic column. In their case, they observed a considerable effect of peak tailing due to mass transfer resistances, and showed the inability of a 1D model to reproduce their results. [START_REF] Augier | Application of computational fluid dynamics to fixed bed adsorption calculations: effect of hydrodynamics at laboratory and industrial scale[END_REF] showed that hydrodynamics inside industrial adsorbers can strongly impair the efficiency of adsorption through the use of CFD. [START_REF] Zheng | Two-dimensional modeling of the transport phenomena in the adsorber during pressure swing adsorption process[END_REF] and [START_REF] Liu | Numerical study of flow maldistribution and depressurization strategies in a small-scale axial adsorber[END_REF] developed two-dimensional models that can operate through cyclic simulations to study a Pressure Swing Adsorption (PSA) process. These models allowed to perceive the process hydrodynamics and to propose better adsorber geometries. They observed concentration maldistributions that would not be visible using a one-dimensional model. This work focuses on a SMB process for the p-xylene (PX) purification. The PX is used in the production of monomers of the polyethylene terephthalate (PET), and is obtained in a mixture with its isomers, o-xylene (OX), m-xylene (MX) and ethylbenzene (EB). Its purification is performed using a set of superposed fixed adsorbers, ranging from 15 to 24 for example, and the simulated countercurrent between the adsorbent and the eluent is achieved by the periodic commutation of the inlet and outlet ports throughout the column that comprises the adsorbers. The SMB contains several hydrodynamic defects: the pipes of the distribution and withdrawal system are placed inside the packing, beams are used to support the weight of the beds, there are empty volumes between consecutive beds and the injection is not ideal, as shown in figure 1(see Pavone and Hotier (2010) and [START_REF] Minceva | Modeling and simulation of a simulated moving bed for the separation of p-xylene[END_REF] for a more detailed description).

The plates of the distribution system are divided into panels that withdraw the flow from the bottom of the adsorbers and mix it with an injected feed from the distribution network to send it to the packed bed below (see for example the patent US20110303602 A1, Augier and Hotier, 2011).

Since its conception, numerous major improvements have been made in the field of the SMB modelling, regarding the mass transfer modelling [START_REF] Mendes | Extension of the linear driving force-dusty gas model approximation to include surface or micropore diffusion[END_REF], CSS prediction (Pais et al., 1997a,b), SMB optimal design [START_REF] Storti | Design of optimal operating conditions of simulated moving bed adsorptive separation units[END_REF]Mazzotti et al., 1997a,b) and column dead volume modelling [START_REF] Migliorini | Simulated moving-bed units with extra-column dead volume[END_REF][START_REF] Minceva | Influence of the transfer line dead volume on the performance of an industrial scale simulated moving bed for p-xylene separation[END_REF][START_REF] Ruthven | Effect of dead volumes on the performance of an industrial-scale simulated moving-bed Parex unit for p-xylene purification[END_REF]). Yet, plug flow models (axially dispersion or not) have been generally used to model the flow adsorption columns and their suitability has not being questioned. [START_REF] Augier | Application of computational fluid dynamics to fixed bed adsorption calculations: effect of hydrodynamics at laboratory and industrial scale[END_REF] studied the impact of a round obstacle within a packed bed of adsorbent on the breakthrough curves. By measuring the resulting variance [START_REF] Augier | Application of computational fluid dynamics to fixed bed adsorption calculations: effect of hydrodynamics at laboratory and industrial scale[END_REF] were able to adjust a DPF model to their case. When using most of the adsorbents (as it will be shown in this work) and after adjusting the axial dispersion coefficient, the DPF model is sufficiently accurate for predicting the hydrodynamic heterogeneities mentioned above. However, if other adsorbents with better mass transfer performances are considered, the accuracy of the DPF model is unknown and has to be investigated.

In order to study the SMB with a proper hydrodynamic description, ideally one should perform CFD simulations coupled with adsorption. However, CFD simulations of multicomponent adsorption in a single adsorber bed can take up to 72 hours to converge (8 cores 2.7 GHz, 16GB RAM) and the SMB comprises a high number of beds. This makes CFD an inappropriate tool for engineers to study multibed cyclic processes. Hence, it is interesting to develop a 1D hydrodynamic model simple enough to converge in seconds but capable of giving a fine description of the adsorbers hydrodynamics. In a previous paper [START_REF] Gomes | Hydrodynamic modelling of complex fixed bed geometries in simulated moving bed adsorption processes[END_REF] developed a 1D model using the CFD methodology developed by [START_REF] Liu | Spatial distributions of mean age and higher moments in steady continuous flows[END_REF] and [START_REF] Liu | Age distribution and the degree of mixing in continuous flow stirred tank reactors[END_REF] based on the work of [START_REF] Danckwerts | The effect of incomplete mixing on homogeneous reactions[END_REF] and [START_REF] Zwietering | The degree of mixing in continuous flow systems[END_REF]. This model has shown a good agreement with the CFD when comparing the inert tracer concentration profiles.

In this work, the accuracy of this 1D model (here called 2MENCSTR model, see abbreviations) coupled with adsorption is evaluated relatively to the DPF model. To do this a chromatographic separation is simulated in a single adsorber. The results obtained with the DPF and the 2MENCSTR model are compared to those given by CFD. Then the DPF and 2MENCSTR models are integrated in a cyclic solver in order to evaluate the importance of a detailed hydrodynamic description when modelling SMB processes.

MODELLING

Adsorption column geometry

In both CFD and 1D models, the following hydrodynamic heterogeneities are taken into account: the packing of adsorbent, a pipe of the distribution/withdrawal system, a beam, and two free flow volumes, one above the packed bed that distributes the flow and one below that withdraws the flow. The dimensions of the geometry used in the CFD model are presented in figure 2.

In this work the flow is considered to be essentially two-dimensional which is a reasonable approximation as shown by [START_REF] Augier | Application of computational fluid dynamics to fixed bed adsorption calculations: effect of hydrodynamics at laboratory and industrial scale[END_REF], allowing much shorter simulation times than if three dimensions were considered. Nevertheless, the methodology proposed can also be used in a threedimensional geometry.

The studied geometry comprises an adsorbent packed bed with 1.2 m height and 0.5 m width, where the half of a pipe (0.151 m radius) and half of a beam (0.268 m edge) are placed. The pipe and beam dimensions were chosen so that each one would occupy around 6% of the packed bed volume. Two free flow volumes of 0.015 m height are placed above and below the packed bed. The inlet and outlet have a length of 0.025 m. The dashed lines in figure 2 represent the symmetry boundary conditions used for both transport and momentum balance.

It is known that the porosity profiles near the walls are variable, however this variation occurs in a range of few particles that constitute the bed. A preliminary study showed that, since the domain in study contains a massive amount of particles, the results obtained were insensible to whether the porosity is assumed to be constant or variable using the profiles proposed by [START_REF] Giese | Measured and modeled superficial flow profiles in packed beds with liquid flow[END_REF].

2.2

Intra-particle phenomena In the adsorbent packed bed, three scales of porosity are used to describe the mass balance: the interstitial porosity 𝜀𝜀, the macro-porosity 𝜀𝜀 𝑝𝑝 and the micro-porosity (or adsorbed phase). The transition of the molecules between the bulk phase and the adsorbed phase is described by two linear driving forces (LDF) [START_REF] Ruthven | Effect of dead volumes on the performance of an industrial-scale simulated moving-bed Parex unit for p-xylene purification[END_REF]. A given component 𝑖𝑖 present in the bulk phase (the concentration of this molecule is expressed as 𝐶𝐶 𝑖𝑖 ) faces a fluid film mass transfer resistance (𝑘𝑘 1 ) when migrating into the macro-porous liquid phase (𝐶𝐶 𝑚𝑚,𝑖𝑖 ). Then, while migrating from the macro-porous liquid phase to the adsorbed phase (𝑞𝑞 𝑖𝑖 ) a second mass transfer resistance is felt (𝑘𝑘 2 ), function of the radius of the adsorbent crystal radius.

Hence, the mass flux of a component 𝑖𝑖 between the bulk phase and the macro-porous phase is given by:

𝑁𝑁 𝑖𝑖 𝜀𝜀 = -(1 -𝜀𝜀)𝜀𝜀 𝑝𝑝 𝑘𝑘 1 �𝐶𝐶 𝑖𝑖 -𝐶𝐶 𝑚𝑚,𝑖𝑖 � (1)
The mass balance of a component 𝑖𝑖 in the macro-porous liquid phase is given by:

𝜀𝜀 𝑝𝑝 𝜕𝜕𝐶𝐶 𝑚𝑚,𝑖𝑖 𝜕𝜕𝜕𝜕 = 𝜀𝜀 𝑝𝑝 𝑘𝑘 1 �𝐶𝐶 𝑖𝑖 -𝐶𝐶 𝑚𝑚,𝑖𝑖 � -�1 -𝜀𝜀 𝑝𝑝 �𝑘𝑘 2 (𝑞𝑞 𝑖𝑖 * -𝑞𝑞 𝑖𝑖 )𝜌𝜌 𝑝𝑝 (2)
for which 𝜌𝜌 𝑝𝑝 is the apparent adsorbent mass density, 𝑞𝑞 𝑖𝑖 * is the adsorbed concentration of 𝑖𝑖 at equilibrium with the macro-porous liquid phase and 𝑞𝑞 𝑖𝑖 is the average adsorbed concentration. The instant average adsorbed concentration of 𝑖𝑖 can be obtained using the adsorbed phase mass balance:

𝜕𝜕𝑞𝑞 𝑖𝑖 𝜕𝜕𝜕𝜕 = 𝑘𝑘 2 (𝑞𝑞 𝑖𝑖 * -𝑞𝑞 𝑖𝑖 ) (3)
The adsorbed concentration at equilibrium with the macro-porous liquid phase can be calculated by the generalized Langmuir isotherm:

𝑞𝑞 𝑖𝑖 * 𝑞𝑞 𝑚𝑚,𝑖𝑖 = 𝑏𝑏 𝑖𝑖 𝐶𝐶 𝑚𝑚,𝑖𝑖 1 + ∑ 𝑏𝑏 𝑗𝑗 𝐶𝐶 𝑚𝑚,𝑗𝑗 𝑁𝑁𝑁𝑁 𝑗𝑗=1 (4)
for which 𝑞𝑞 𝑚𝑚,𝑖𝑖 is the adsorbent mass capacity for the component 𝑖𝑖, 𝑏𝑏 𝑖𝑖 is the adsorbent coefficient of affinity for the component 𝑖𝑖 and 𝑁𝑁𝐶𝐶 is the number of components in equilibrium.

CFD Model

Even though adsorption is often an exothermic phenomenon, in this study the adsorbent is permanently saturated and the molecules show similar values of adsorption capacity and enthalpy. Thus, it is reasonable to consider it as isothermal. With this, the mass and momentum balance are calculated but the energy balance is not.

In a previous paper [START_REF] Gomes | Hydrodynamic modelling of complex fixed bed geometries in simulated moving bed adsorption processes[END_REF] found a good agreement between the concentration profiles of an inert tracer obtained experimentally and numerically, giving a fair degree of thrust in the current CFD model (COMSOL 4.3b software).

Since the geometry in study comprises a packed bed between two free flow volumes, the momentum and mass balance must be solved differently leading to two distinct sets of equations.

Given the fact that at the SMB operating conditions (180 °C and 9 bar, [START_REF] Minceva | Understanding and revamping of industrial scale SMB units for p-xylene separation[END_REF] the xylenes are in liquid phase thus incompressible and that the inlet and outlet volumetric flows are constant the flow is assumed to be steady in both porous and free flow.

Free flow

In the free flow volumes placed above and below the packed bed the flow is described using the Reynolds averaged Navier-Stokes equations (RANS) for an incompressible Newtonian fluid and steady flow. The turbulent dynamic viscosity 𝜇𝜇 𝑇𝑇 is modelled by the standard k-ε turbulence model [START_REF] Jones | The prediction of laminarization with a two-equation model of turbulence[END_REF].

With such macroscopic description of the flow, the mixing of the molecules due to eddy motion is not simulated. This phenomenon is then modelled by imposing a constant turbulent Schmidt number:

𝑆𝑆𝑆𝑆 𝑇𝑇 = 𝜇𝜇 𝑇𝑇 𝜌𝜌𝐷𝐷 𝑇𝑇 ( 5
)
𝜇𝜇 𝑇𝑇 is the turbulent dynamic viscosity calculated by the k-ε turbulence model, 𝜌𝜌 is the fluid mass density and 𝐷𝐷 𝑇𝑇 is the turbulent diffusivity coefficient. This diffusivity coefficient must be added to the molecular diffusion 𝐷𝐷 𝑀𝑀 on the diffusive term of the mass balance equation, in the free flow regions. Based on the average values found in the bibliography [START_REF] Koeltzsch | The height dependence of the turbulent Schmidt number within the boundary layer[END_REF][START_REF] Flesch | Turbulent Schmidt number from a tracer experiment[END_REF][START_REF] Combest | On the gradient diffusion hypothesis and passive scalar transport in turbulent flows[END_REF] the 𝑆𝑆𝑆𝑆 𝑇𝑇 is assumed to be 0.7.

In the free flow volumes a component 𝑖𝑖 is transported by diffusion and convection:

𝜕𝜕𝐶𝐶 𝑖𝑖 𝜕𝜕𝜕𝜕 = ∇ • �(𝐷𝐷 𝑀𝑀 + 𝐷𝐷 𝑇𝑇 ) • ∇𝐶𝐶 𝑖𝑖 -𝑢𝑢 �⃗𝐶𝐶 𝑖𝑖 � (6)
The global dispersion coefficient results from the sum of the molecular diffusion and the turbulent diffusivity, both isotropic.

Porous media

The CFD model describes hydrodynamics at a scale much larger than the bed particles which is reasonable given the high ratio between the bed dimensions and the particles diameter. The interaction between the fluid and the particles is not simulated, but modelled by the addition of a diffusive term in the momentum and mass balance. Then, in order to simulate the flow in the porous media a modification of the laminar Navier-Stokes equations is used, called the Brinkman-Forchheimer equations [START_REF] Vafai | Boundary and inertia effects on convective mass transfer in porous media[END_REF]:

𝜌𝜌 𝜀𝜀 �(𝑢𝑢 �⃗ • ∇) 𝑢𝑢 �⃗ 𝜀𝜀 � = ∇ • �-𝑝𝑝 + 𝜇𝜇 𝑒𝑒 𝜀𝜀 (∇𝑢𝑢 �⃗ + (∇𝑢𝑢 �⃗) 𝑇𝑇 ) - 2𝜇𝜇 3𝜀𝜀 (∇ • 𝑢𝑢 �⃗)� -�𝜌𝜌 • 𝛽𝛽 • |𝑢𝑢 �⃗| + 𝜇𝜇 𝐾𝐾 � 𝑢𝑢 �⃗ (7) 𝜌𝜌 • ∇𝑢𝑢 �⃗ = 0 (8)
𝑢𝑢 �⃗ is the superficial fluid velocity, 𝜀𝜀 the bed porosity, p the pressure, 𝜇𝜇 𝑒𝑒 the effective dynamic viscosity, 𝜇𝜇 the dynamic viscosity, 𝛽𝛽 is the non-Darcy coefficient and 𝐾𝐾 the bed permeability. The validity of these equations has been discussed by several authors [START_REF] Nield | The limitations of the Brinkman-Forchheimer equation in modeling flow in a saturated porous medium and at an interface[END_REF][START_REF] Vafai | On the limitations of the Brinkman-Forchheimer-extended Darcy equation[END_REF][START_REF] Chan | Numerical study of forced flow in a back-step channel through porous layer[END_REF][START_REF] Benyahia | On the modeling of flow in packed bed systems[END_REF][START_REF] Augier | Application of computational fluid dynamics to fixed bed adsorption calculations: effect of hydrodynamics at laboratory and industrial scale[END_REF][START_REF] Gomes | Hydrodynamic modelling of complex fixed bed geometries in simulated moving bed adsorption processes[END_REF]. The value of 𝜇𝜇 𝑒𝑒 has been a subject of debate, but similarly to the previous paper [START_REF] Gomes | Hydrodynamic modelling of complex fixed bed geometries in simulated moving bed adsorption processes[END_REF], where the validity of the current CFD model is evaluated, its value is assumed to be equal to 𝜇𝜇.

𝛽𝛽 and 𝐾𝐾 can calculated, for a bed of spheres by the Ergun's law (1952):

𝛽𝛽 = 1.75 1 -𝜀𝜀 𝜀𝜀 3 𝑑𝑑 𝑝𝑝 (9) 𝐾𝐾 = 𝜀𝜀 3 𝑑𝑑 𝑝𝑝 2 150(1 -𝜀𝜀) 2 (10)
As aforementioned, a diffusive term must be added to the mass balance equation in order to model the dispersion caused by the packing. This is done by adding the mechanical dispersion coefficient 𝐷𝐷 𝑃𝑃 to the molecular diffusion coefficient, in the diffusive term of the mass balance equation of the bulk phase:

𝜀𝜀 𝜕𝜕𝐶𝐶 𝑖𝑖 𝜕𝜕𝜕𝜕 + ∇ • (𝑢𝑢 �⃗𝐶𝐶 𝑖𝑖 -𝜀𝜀(𝐷𝐷 𝑀𝑀 + 𝐷𝐷 𝑃𝑃 ) • ∇𝐶𝐶 𝑖𝑖 ) = 𝑁𝑁 𝑖𝑖 𝜀𝜀 (11)
The mass flux of a component 𝑖𝑖 between the bulk phase and the macro-porous phase, 𝑁𝑁 𝑖𝑖 𝜀𝜀 , is given by Eq. ( 1).

Contrarily to the molecular diffusion (𝐷𝐷 𝑀𝑀 ) and turbulent diffusivity (𝐷𝐷 𝑇𝑇 ) coefficients, the mechanical dispersion coefficient (𝐷𝐷 𝑃𝑃 ) is anisotropic and described by a diagonal tensor containing the transversal and axial dispersion coefficients. These are calculated by imposing two constant Peclet numbers based on the particle diameter, one in the transversal direction (𝑃𝑃𝑃𝑃 𝑡𝑡𝑡𝑡 ) and the other in the axial direction (𝑃𝑃𝑃𝑃 𝑎𝑎𝑎𝑎 ), relatively to the flow direction. The values of 𝑃𝑃𝑃𝑃 𝑡𝑡𝑡𝑡 and 𝑃𝑃𝑃𝑃 𝑎𝑎𝑎𝑎 are assumed to be 11 and 2 respectively, the asymptotic values for liquid flows with a particle Reynolds number higher than 20 [START_REF] Delgado | A critical review of dispersion in packed beds[END_REF].

To couple the two free flow volumes with the porous media, the boundary conditions in use force the pressure, the horizontal and the vertical components of the fluid velocity and the bulk phase concentration of every component to be the same on both sides of the interface, for each mesh element.

2.4

Degree of mixing and transport of the fluid age distribution moments In 1952, Danckwerts showed that two reactors that share the same Residence Time Distribution (RTD) can produce different conversions, if the chemical reaction follows a second order kinetics. To explain this, [START_REF] Danckwerts | The effect of incomplete mixing on homogeneous reactions[END_REF] and [START_REF] Zwietering | The degree of mixing in continuous flow systems[END_REF] formulated a criterion, the degree of mixing (𝐽𝐽), that quantifies the mixing earliness of a vessel. However, for about 50 years it was not possible to calculate its value, except for the asymptotic values: 0 for an ideal Continuous Stirred-tank Reactor (CSTR) and 1 for a Plug Flow Reactor (PFR). Recently, [START_REF] Liu | Age distribution and the degree of mixing in continuous flow stirred tank reactors[END_REF] showed that the degree of mixing could be computed using the volume averages of the first and second raw moments of the fluid age distribution (𝑚𝑚 1 and 𝑚𝑚 2 ) whose spatial distribution can be obtained by coupling a steady moment balance equation with the flow field simulated by CFD [START_REF] Liu | Spatial distributions of mean age and higher moments in steady continuous flows[END_REF]. Such methodology can only be applied to steady incompressible flows.

Adsorption is a non-linear exchange phenomena in which 𝑞𝑞 𝑖𝑖 * varies non-linearly with 𝐶𝐶 𝑚𝑚,𝑖𝑖 . Therefore, when proposing a 1D model capable of reproducing the coupling between hydrodynamics and adsorption, it is of the highest interest to preserve the CFD model degree of mixing. The degree of mixing obtained for the adsorber geometry represented in figure 2 is close to 1 (𝐽𝐽=0.987), indicating a behaviour close to a completely segregated system, despite the hydrodynamic heterogeneities [START_REF] Gomes | Hydrodynamic modelling of complex fixed bed geometries in simulated moving bed adsorption processes[END_REF].

1D Model

As long as the degree of mixing is close to the one of a PFR it is reasonable to assume that no transversal mixing occurs in the adsorber. For this reason, each trajectory of the fluid is modelled by three consecutive cascades of 𝑁𝑁 CSTR models (hereby called NCSTR models), one for each zone (above the packed bed, below the bed and the bed itself), as shown in figure 3 a). Finally, when all the fluid trajectories are modelled in this way, the CFD model can be condensed to a set of parallel systems comprising three consecutive NCSTR models, as shown in figure 3 b [START_REF] Gomes | Hydrodynamic modelling of complex fixed bed geometries in simulated moving bed adsorption processes[END_REF]. [START_REF] Gomes | Hydrodynamic modelling of complex fixed bed geometries in simulated moving bed adsorption processes[END_REF] observed that the dispersive effect of the zone above the bed (zone A) is induced by the heterogeneity of 𝑚𝑚 1 in the inlet of the porous media (the fluid elements take different times to flow from the inlet to the entrance of the packed bed) rather than by local dispersion. This fact allows to model both free flow volumes downstream of the bed, simplifying the model. One major assumption that is made in order to develop this model is that the horizontal velocity profiles are quite flat, i.e. no major velocity variations are found in the x-axis.

A given streamline, modeled by a cascade of N CSTR has a first portion identical to the cascade of M CSTR corresponding to a smaller streamline. Therefore both streamlines can be modeled with the aid of a single CSTR cascade comprising two different parts. A first part consists in a cascade of M CSTR in which the flow passing through is equal to the sum of the flows of both streamlines. At the end of this first part, the flow corresponding to the shortest streamline is withdrawn, and the flow corresponding to the longest streamline flows through a cascade of (N-M) CSTR. By doing this individually for each streamline, it is possible to represent the porous zone with a Multiple Exit NCSTR model and the free flow zones with a Multiple Entrance NCSTR model, henceforth denominated 2MENCSTR model (as for Multiple Exits NCSTR model connected to a Multiple Entrances NCSTR model), as shown in figure 3 c.

SMB modelling

In order to model a SMB process the switching of the inlet and outlet ports is accounted for, as shown in figure 4 a, using a time-dependent solver. The solid arrows represent the streams at an instant 𝜕𝜕 while the dashed arrows represent the streams after the stream switching, at 𝜕𝜕+∆𝜕𝜕 𝑠𝑠𝑠𝑠 , where ∆𝜕𝜕 𝑠𝑠𝑠𝑠 is the switching time.

Extract node

𝑄𝑄 𝑍𝑍2 = 𝑄𝑄 𝑍𝑍1 -𝑄𝑄 𝑋𝑋 (14) 𝐶𝐶 𝑖𝑖,𝑋𝑋 = 𝐶𝐶 𝑖𝑖,𝑗𝑗+1 𝑖𝑖𝑖𝑖 = 𝐶𝐶 𝑖𝑖,𝑗𝑗 𝑜𝑜𝑜𝑜𝑡𝑡 (15)
Feed node

𝑄𝑄 𝑍𝑍3 = 𝑄𝑄 𝑍𝑍2 + 𝑄𝑄 𝐹𝐹 (16) 𝑄𝑄 𝑍𝑍3 𝐶𝐶 𝑖𝑖,𝑗𝑗+1 𝑖𝑖𝑖𝑖 = 𝑄𝑄 𝑍𝑍2 𝐶𝐶 𝑖𝑖,𝑗𝑗 𝑜𝑜𝑜𝑜𝑡𝑡 + 𝑄𝑄 𝐹𝐹 𝐶𝐶 𝑖𝑖,𝐹𝐹 (17)
Raffinate node

𝑄𝑄 𝑍𝑍4 = 𝑄𝑄 𝑍𝑍3 -𝑄𝑄 𝑅𝑅 (18) 𝐶𝐶 𝑖𝑖,𝑅𝑅 = 𝐶𝐶 𝑖𝑖,𝑗𝑗+1 𝑖𝑖𝑖𝑖 = 𝐶𝐶 𝑖𝑖,𝑗𝑗 𝑜𝑜𝑜𝑜𝑡𝑡 (19)
where 𝑄𝑄 𝑍𝑍𝑖𝑖 is the volumetric flow rate of zone 𝑖𝑖, the index 𝐷𝐷 stands for desorbent, 𝑋𝑋 for extract, 𝐹𝐹 for feed and 𝑅𝑅 for raffinate. By employing these equations and periodically commuting the inlet and outlet nodes through the SMB column during several cycles it is possible to achieve a cyclic steady state (CSS).

The equation systems reported above are implemented in a time-dependent solver in order to model the separation of a binary mixture of p-xylene (PX) and m-xylene (MX), using p-diethylbenzene (PDEB) as desorbent. Then, after the CSS is reached, the concentration profiles and performances obtained by the DPF-based SMB model and the 2MENCSTR-based SMB model are compared, for the same geometry and operating conditions.

Computation

The CFD simulations were performed with COMSOL 4.3b software. The flow field was solved using a quadratic discretization for the pressure and linear discretization for the velocity. The mass balance equations were solved using a quadratic discretization and a second order BDF for the time stepping.

The 1D models (DPF and 2MENCSTR) computing code has been developed in Fortran (including the SMB solver). The ordinary differential equations system obtained by spatial discretization by cascades of CSTRs is solved using the DDASPK subroutine based on the Petzold-Gear BDF method.

Simulations parameters

The parameters used in this study are shown in Table 1. The equilibrium values 𝑞𝑞 𝑖𝑖 and 𝑏𝑏 𝑖𝑖 [START_REF] Minceva | Understanding and revamping of industrial scale SMB units for p-xylene separation[END_REF] and the LDF constants 𝑘𝑘 1 and 𝑘𝑘 2 [START_REF] Bergeot | Intensification of paraxylene production using a simulated moving bed reactor[END_REF] were found in the bibliography. For the chromatographic simulations, the inlet velocity was adjusted with the aim of setting the residence time inside the porous bed to 27.4 s, which leads to an average superficial velocity of 1.4 cm/s, close to what is used in industrial units. It must be reminded that the dimensionless numbers Schmidt and Peclet were not fitted but obtained from the bibliography. Two different SMB configurations were considered in this work, comprising 6 and 12 adsorption columns. The SMB operating conditions used are shown in Table 2.

RESULTS AND DISCUSSION

DTS and model parameterization

The characteristic parameters of the resulting 2MENCSTR model, that represents the geometry shown in Fig. 2, are shown in figure 5.

In figure 5, the first column shows then cumulative number of CSTR in the first cascade of the 2MENCSTR model, i.e. the porous one. The number of CSTR is counted until a side exit, i.e. the first side exit comes out from the 1846 th CSTR, while the second side exit comes out from the 1852 th CSTR, and so on. The second column shows the total volume of each CSTR from a given excerpt of the cascade, including the volume of the packing and the void fraction. All CSTR between 1 and 1845 have the same volume that the 1846 th CSTR.

Similarly, the first and second rows show the number and volumes of the CSTR in the second cascade. Within the matrix, the flows exchanged by the two cascades are represented as fractions of the flow rate at the inlet of the first CSTR of the bed.

As mentioned above, the single arbitrary parameter of this model is the number of interconnections between the two cascades. In this study, the number of interconnections was set to 20, since it is enough to give a good temporal discretization of the adsorption column RTD while keeping the model simple.

Concerning the DPF model, the axial dispersion coefficient is calculated through the RTD moments (𝑃𝑃𝑃𝑃 2 ⁄ = 𝜇𝜇 1 2 𝜎𝜎 2 ⁄ with a Peclet number involving the bed length and not the particle diameter). For simplification, the equivalent DPF model is assimilated to a cascade of 73 CSTR (corresponding to a Peclet number of 146), with a higher interstitial porosity (33.878%) in order to take into account of the free flow zones above and below the bed.. The inlet flow considered for RTD simulation is 3600 m3/h. At t=0, a Dirac injection modelled by a Gaussian shape (𝑃𝑃 -50(𝑡𝑡-1) 2 where 𝜕𝜕 is time) is introduced. RTD calculated with CFD, DPF and 2MENCSTR models are compared in figure 6. The asymmetry of the CFD curve is well reproduced by the 2MENCSTR. The DPF model generates a symmetric RTD.

3.2

Chromatographic separation In order to compare the 2MENCSTR model with the CFD and simultaneously evaluate the dispersive effect produced by each one of the hydrodynamic heterogeneities, in an initial state the adsorber is filled with p-diethylbenzene (PDEB), the desorbent regularly used in the PX purification. Then, a Dirac injection modelled by a Gaussian shape (𝑃𝑃 -50(𝑡𝑡-1) 2 where 𝜕𝜕 is time) of a mixture comprising 95% of PX and 5% of MX is introduced.

In figure 7, the concentration profiles of PX and MX obtained using the parameters shown in Table 1 in the CFD model, a DPF model and the model proposed, the 2MENCSTR model are shown.

The deviations between the CFD and the DPF models are very slight but the 2MENCSTR fits better the CFD results. In this case the use of the 2MENCSTR model, instead of the DPF model, does not present a significant gain regarding the coupling between hydrodynamics and adsorption. Since the RTD of a DPF model is symmetric, and that the concentration profiles obtained with this model are visibly leaning to the right, it is possible to conclude that in this case most of the peak tailing noticed is due to the slow mass transfer rates, rather than the adsorber design and the consequent flow patterns.

In order to investigate if this SMB-enhancing strategy is worth and also to further test the models in study, in a scenario where the adsorbent in use has a better mass transfer capacity, the LDF constants presented in Table 1 are multiplied by a factor of 10. The resulting concentration profiles are shown in figure 8.

While the DPF model is unable to reproduce the asymmetric concentration profiles of the CFD model, the 2MENCSTR shows a good agreement with the CFD, yet showing some deviations when the concentration profiles are compared in a logarithmic scale.

In order to investigate the importance of a detailed hydrodynamic description in function of the mass transfer scenario, the chromatographic separation is simulated with the three models in study, CFD, DPF and 2MENCSTR for five different sets of LDF constants 𝑘𝑘 1 and 𝑘𝑘 2 . The LDF constants used by [START_REF] Bergeot | Intensification of paraxylene production using a simulated moving bed reactor[END_REF] (𝑘𝑘 1 =0.5 s -1 and 𝑘𝑘 2 =0.153 s -1

) were multiplied by 0.5, 1, 5, 10 and 50 and the deviation between the 1D model (either DPF or 2MENCSTR) and the CFD model was measured by equation 20 and reported on figure 9.

𝑑𝑑𝑃𝑃𝑑𝑑𝑖𝑖𝑑𝑑𝜕𝜕𝑖𝑖𝑑𝑑𝑑𝑑 = � ∫ (𝐶𝐶 𝑃𝑃𝑋𝑋 𝑁𝑁𝐹𝐹𝐷𝐷 -𝐶𝐶 𝑃𝑃𝑋𝑋 1𝐷𝐷 ) 2 𝑑𝑑𝜕𝜕 ∞ 0 ∫ 𝐶𝐶 𝑃𝑃𝑋𝑋 𝑁𝑁𝐹𝐹𝐷𝐷 2 𝑑𝑑𝜕𝜕 ∞ 0 (20)
Figure 9 shows that when an adsorbent with mass transfer performances equal or lower than the one used by [START_REF] Bergeot | Intensification of paraxylene production using a simulated moving bed reactor[END_REF] is used, the DPF model is capable of predicting the concentration profiles with a reasonable error. The DPF model shows a deviation from the CFD concentration profiles of 7%, slightly higher than the 2MENCSTR model with 5% for k LDF multipliers of 0.5 and 1. A fine description of hydrodynamics is not mandatory. However, if the adsorbent used is capable of inducing mass transfer rates faster than those indicated before (k LDF multiplier of 5, 10 and 50), the deviations between the DPF and the CFD model are more important than the deviations between the 2MENCSTR and the CFD (that also increase with the k LDF multiplier). For instance, for a k LDF multiplier of 10, the DPF model shows a deviation of 27% from the CFD concentration profiles while the 2MENCSTR deviation is 3 times lower. In the case of fast mass transfer rates, it is interesting to use a more detailed 1D hydrodynamic model in order to minimize the discrepancies between a CFD model (2D in this case but also applicable to 3D geometries) and a simpler and faster model.

It is then acceptable that the 2MENCSTR model has been validated for SMB studies by comparing its results with those obtained with a CFD model of a single adsorption column. The 2MENCSTR model is then integrated in a cyclic solver and compared to the DPF model to evaluate the usefulness of a detailed hydrodynamic description in SMB modelling.

SMB results

The PX and MX profiles obtained for a SMB unit comprising 6 adsorption columns are shown in the form of mass fraction for 𝑘𝑘 1 and 𝑘𝑘 2 equal to 5 and 1.53 s, respectively, at linear and logarithmic scales in figure 10.

In order to quantify the SMB performance, the purity and the recovery of PX in the extract stream are measured:

𝑃𝑃𝑃𝑃𝑅𝑅 𝑋𝑋 = 〈𝐶𝐶 𝑃𝑃𝑋𝑋,𝑋𝑋 〉 〈𝐶𝐶 𝑃𝑃𝑋𝑋,𝑋𝑋 〉 + 〈𝐶𝐶 𝑀𝑀𝑋𝑋,𝑋𝑋 〉 (21)
𝑅𝑅𝑅𝑅𝐶𝐶 𝑋𝑋 = 𝑄𝑄 𝑋𝑋 〈𝐶𝐶 𝑃𝑃𝑋𝑋,𝑋𝑋 〉 𝑄𝑄 𝐹𝐹 𝐶𝐶 𝑃𝑃𝑋𝑋,𝐹𝐹 (22)
〈𝐶𝐶 𝑃𝑃𝑋𝑋,𝑋𝑋 〉 and 〈𝐶𝐶 𝑀𝑀𝑋𝑋,𝑋𝑋 〉 denote that the concentration of PX and MX are averaged on the extract stream over a switching time (after the CSS is reached), due to the SMB time-variant concentration profiles.

When comparing the concentration profiles obtained with the 2MENCSTR and DPF models, slight deviations are visible, both at linear and logarithmic scales. In order to quantify these deviations the 2MENCSTR model is taken as a reference and the deviations between the extract purity and recovery provided by both models is quantified. These deviations are more penalizing at high purities and recoveries, i.e. an absolute deviation of 1% is more prejudicial for a purity of 99% than for a purity of 90%. Thus, the deviations between both models are calculated in terms of relative deviations from 100%, and reported in table 3:

-considering the deviation of impurities fraction in the extract:

𝜖𝜖 𝑃𝑃𝑃𝑃𝑅𝑅 𝑋𝑋 = (1 -𝑃𝑃𝑃𝑃𝑅𝑅 𝑋𝑋 𝐷𝐷𝑃𝑃𝐹𝐹 ) -(1 -𝑃𝑃𝑃𝑃𝑅𝑅 𝑋𝑋 2𝑀𝑀𝑀𝑀𝑁𝑁𝑁𝑁𝑀𝑀𝑇𝑇𝑅𝑅 ) (1 -𝑃𝑃𝑃𝑃𝑅𝑅 𝑋𝑋 2𝑀𝑀𝑀𝑀𝑁𝑁𝑁𝑁𝑀𝑀𝑇𝑇𝑅𝑅 ) (23)
-considering the deviation of PX fraction in the raffinate:

𝜖𝜖 𝑅𝑅𝑀𝑀𝑁𝑁 𝑋𝑋 = (1 -𝑅𝑅𝑅𝑅𝐶𝐶 𝑋𝑋 𝐷𝐷𝑃𝑃𝐹𝐹 ) -(1 -𝑅𝑅𝑅𝑅𝐶𝐶 𝑋𝑋 2𝑀𝑀𝑀𝑀𝑁𝑁𝑁𝑁𝑀𝑀𝑇𝑇𝑅𝑅 ) (1 -𝑅𝑅𝑅𝑅𝐶𝐶 𝑋𝑋 2𝑀𝑀𝑀𝑀𝑁𝑁𝑁𝑁𝑀𝑀𝑇𝑇𝑅𝑅 ) (24)
By looking at Table 3, it is possible to observe that the DPF model strongly underestimates the PX purity and recovery when compared to the 2MENCSTR model that is used as a reference. The deviation of residual impurities in the extract line is over 80%. This signifies that using the DPF model to design a SMB process may lead to a loss of productivity and/or yield.

In order to further investigate the importance of a detailed hydrodynamic description in SMB modelling, the same study was conducted for a 12 column SMB process. The concentration profiles obtained are presented in figure 11.

Even though in this case, the concentration profiles provided by both models appear to provide a better agreement than in the case of the SMB with 6 adsorption columns, a major profile deviation is found near the extract (bed 2) stream at both linear and logarithmic scales. The resulting performance parameters are shown in table 4.

The DPF model underestimates the extract purity, as for the 6 column SMB, with an absolute deviation of 0.21% and a deviation of residual impurities of 210.00% , while the extract recovery deviations are of 0.43% (in absolute) and 70.49% (relative loss of PX in the raffinate line). When these values are compared to those obtained for a 6 bed SMB, the absolute deviations decrease, as the SMB leans towards a "perfect operation" and both extract purity and recovery approach 100%. However the deviations calculated by Eqs. 23 and 24 are more important for the 12 bed SMB, meaning that a detailed hydrodynamic description is more important as the SMB approaches high purities and recoveries. [START_REF] Gomes | Hydrodynamic modelling of complex fixed bed geometries in simulated moving bed adsorption processes[END_REF] proposed a method to develop 1D models based on the transport of the moments of the fluid age distribution by CFD [START_REF] Liu | Age distribution and the degree of mixing in continuous flow stirred tank reactors[END_REF][START_REF] Liu | Spatial distributions of mean age and higher moments in steady continuous flows[END_REF][START_REF] Liu | Age distribution and the degree of mixing in continuous flow stirred tank reactors[END_REF]. In this work, this method is used and the 1D model deduced can simulate an adsorption process more accurately than the traditionally used DPF model, when coupled with adsorption. The accuracy gain when using the 1D model proposed, the 2MENCSTR model, instead of the DPF model is more evident when the mass transfer rates are fast. The 2MENCSTR can converge in few minutes, contrasting with the CFD simulations of multicomponent adsorption coupled with hydrodynamics that can take several days to converge at present. This facilitates the parametric study and optimization of the SMB or other industrial applications involving adsorbers.

CONCLUSION

The adsorber geometry reveals itself to be a minor issue for slow mass transfer rates, yet not negligible. However, if the development of an adsorbent with better mass transfer performances is considered the adsorber geometry becomes preponderant. The addition of free flow volumes and obstacles to the packed bed can severely decrease the separation efficiency of the adsorber.

Then, the 2MENCSTR model was used as a reference in a cyclic solver and compared to the DPF model to investigate the importance of a good hydrodynamic description when studying SMB processes. The DPF model is able to provide a reasonable agreement regarding the 2MENCSTR SMB concentration profiles but the performance parameters are visibly different from the ones obtained with the 2MENCSTR model, reaching important extract purity and recovery deviations when residual impurity concentrations are compared.

In this work 1D models were deduced for two-dimensional adsorber geometries, although the same method can be used for 3D geometries and for diverse applications of separation processes or reactors.

NOTATIONS

𝑏𝑏 𝑖𝑖

adsorbent coefficient of affinity for the component 𝑖𝑖 (m Table 3 -PX purity and recovery obtained with the DPF and 2MENCSTR models for a SMB unit with 6 adsorption columns, for k 1 and k 2 equal to 5 and 1.53 s, respectively.

Table 4 -PX purity and recovery obtained with the DPF and 2MENCSTR models for a SMB unit with 12 adsorption columns, for k 1 and k 2 equal to 5 and 1.53 s, respectively.

Figure Captions

Figure 1 -Schematic view of an example of SMB process [START_REF] Gomes | Hydrodynamic modelling of complex fixed bed geometries in simulated moving bed adsorption processes[END_REF].

  Then, two different models are used to model the SMB adsorption columns: the 2MENCSTR model and the DPF model. The 2MENCSTR model ( figure 4 b) comprises a cascade of CSTR models with multiple side exits that are connected to a second CSTR cascade with multiple side entrances. The first cascade is used to model the porous bed, thus the CSTR are filled with adsorbent, while the second cascade is used to model the free flow regions above and below the bed, thus the CSTR are "empty". The parameters regarding the CSTR sizes, the flow of each connection and connection arrangement are function of the spatial distributions of the moments of the fluid age distribution transported by CFD. The single arbitrary parameter is the number of interconnections. DPF models are generally used for the modelling of adsorption columns in SMB studies. In this work the Peclet numbers of the DPF models are adjusted to match the variance of the associated 2MENCSTR model. It is assured that for each SMB configuration and geometry in study the 2MENCSTR and the DPF models have the same residence time, variance and mass of adsorbent. Whether in the use of the 2MENCSTR or DPF model the following equations are used to perform the mass balance of a component 𝑖𝑖 in the inlet and outlet nodes between two consecutive adsorption columns 𝑗𝑗 and 𝑗𝑗+1: Desorbent node 𝑄𝑄 𝑍𝑍1 = 𝑄𝑄 𝑍𝑍4 + 𝑄𝑄 𝐷𝐷 (12) 𝑄𝑄 𝑍𝑍1 𝐶𝐶 𝑖𝑖,𝑗𝑗+1 𝑖𝑖𝑖𝑖 = 𝑄𝑄 𝑍𝑍4 𝐶𝐶 𝑖𝑖,𝑗𝑗 𝑜𝑜𝑜𝑜𝑡𝑡 + 𝑄𝑄 𝐷𝐷 𝐶𝐶 𝑖𝑖,𝐷𝐷

Figure 2 -

 2 Figure 2 -Geometry and dimensions of the adsorption bed.

Figure 3 -

 3 Figure 3 -Scheme of the modelling approach proposed (a), resulting system of parallel NCSTR models with the distinction the dissociation of the 3 zones (b) and resulting 2MENCSTR model (c).

Figure 4 -

 4 Figure 4 -(a) Schematic diagram of a SMB process for the separation of a binary mixture, (b) 2MENCSTR model and (c) DPF model, both used to model each adsorption column.

Figure 5 -

 5 Figure 5 -Matrix with the fractions of flow exchanged between the two CSTR cascades of the 2MENCSTR model, number (N) and volume (V) of the CSTR.

Figure 6 -

 6 Figure 6 -RTD comparison between the different models of similar variances σ².

Figure 7 -

 7 Figure 7 -PX and MX concentration profiles obtained with the three studied models; CFD, DPF and 2MENCSTR.

Figure 8 -

 8 Figure 8 -PX and MX concentration profiles obtained with the three studied models; CFD, DPF and 2MENCSTR, for k 1 and k 2 equal to 5 and 1.53 s, respectively.

Figure 9 -

 9 Figure 9-Deviations between the concentration profiles obtained by CFD and obtained with the 1D models, DPF and 2MENCSTR for different mass transfer scenarii.

Figure 10 -

 10 Figure10-Mass fraction profiles obtained with the DPF and 2MENCSTR models for a SMB with 6 adsorption columns, for k 1 and k 2 equal to 5 and 1.53 s, respectively.

Figure 11 -

 11 Figure11-Mass fraction profiles obtained with the DPF and 2MENCSTR models for a SMB with 12 adsorption columns, for k 1 and k 2 equal to 5 and 1.53 s, respectively.

  Figure 6

Table Captions Table 1 -

 Captions1 Simulation parameters.

Table 2 -

 2 SMB operating conditions.