Structure-preserving infinite dimensional model reduction: Application to adsorption processes

Structure-preserving infinite dimensional model reduction Application to Adsorption Processes. A.Baaiu, F.Couenne, L.Lefevre, Y.Le Gorrec † , M.Tayakout June 3, 2008 LAGEP, UMR C.N.R.S. , University Lyon 1, Bat 308G, 43 Bd du 11 novembre 1918, 69622 Villeurbanne Cedex, France. † Email:[email protected] Abstract This paper deals with the model reduction (spatial discretization) of distributed parameter systems. The originality of the proposed approach is to combine geometrical modelling and finite element discretization method to preserve the model structure associated with both mass and energy balances during the reduction. The approach is presented through the example of an adsorption process whose main feature is that it is an heterogeneous system. The methodology is described on the microporous level (phase) and simulations are carried out on the full multilevel adsorption column.

Keywords : port-based modelling, model reduction, infinite dimensional systems, thermodynamic systems

1

Introduction

An adsorption column is a complex heterogeneous system which can be mathematically described by a set of partial differential equations (PDE’s) [21, 22, 34, 35]. Due to heterogeneity and the use of porous medium, the physical system can be considered respect to three different spatially nested levels (also denoted scale in the text). The traditional modelling of such a system over these three scales is based on balance equations. This modelling doesn’t take into account any structure of the constitutive equations. The choice of the manipulated variables can lead to numerical difficulties, particularly during the interconnection of the different scales due to the non consistency of the boundary conditions. Furthermore, usually the model reduction is obtained using spatial discretization strategies that are unconnected with the modelling strategy. An alternative way to model the adsorption column is to use a network approach, similar to the one used for non dissipative systems [24], to stress the intrinsic structure of the physical phenomena occurring into the column. This network approach consists in splitting the model into a set of atomistic phenomena with characteristic energetic behavior. These atomistic elements are joined together using an interconnection structure acting on some power conjugated variables, named port variables. This interconnection structure characterizes the lossless power exchanges within the system and through its boundaries. The fact that the interconnection of the different levels is done using power conjugate variables induces good properties during the interconnection (equality of efforts and continuity of flux). From an energetic point of view, three basic elements are sufficient to characterize the main energetic phenomena occurring within the column: the storage element, the dissipation element and the interconnection element which captures the distributed aspect of the conservation laws. With this representation are stressed the internal energy variation, power flow at the boundary and passivity properties. These properties are fundamental in control theory since they can be used for stability analysis or control purposes [27, 30, 31, 13]. As structural properties of the model are pointed out, a natural extension of this approach concerns the geometric reduction of distributed parameters systems. Finite dimensional approximation (also called reduction) is a fundamental concern for the simulation and control of distributed parameters complex systems. Traditional methods used 1

for the spatial discretization of PDEs are based on the approximation of equations or the approximation of solutions [39]. The first category gather the well known finite differences method and finite volume method while the second one cops with weighted residual methods, finite element methods etc ... One of the difficulties is to develop a method of reduction that preserves some interesting qualitative features of the original model (such as stability, passivity, etc...). Inspired from the aforementioned modelling strategy, we present in this paper a new method of spatial discretization for infinite dimensional systems which applies to many distributed parameters thermodynamic systems. It is based on the preservation of the local structure of the power exchanges and can be assimilated to mixed finite elements method [29]. The adsorption column model is used as an illustrative example. To summarize, the aim of this paper is to present a new spatial discretization method associated with the network modelling of irreversible infinite dimensional thermodynamic systems with an application to adsorption processes. It is an extension of the work presented in [12] in the case of reversible electrical systems. The paper is organized as follows. In Section 2, we briefly recall some basic theoretical features for the port-based modelling of distributed thermodynamic parameters systems. Such structured modelling is based on the fundamental Gibbs equation expressing the thermodynamic properties of a mixture, the balance equations and some closure equations related to the considered phenomena. This section is quite general and gives motivation for the geometrical modeling of these systems. From mathematical point of view, the use of the differential geometry allows to obtain coordinates free models. It makes the notations and structure general and applicable to 1-D, 2-D or 3-D systems. The concrete application to the adsorption case is given in Section 3. In this section we first introduce the modelling of adsorption processes under consideration and present the associated with multilevel geometric model. In this section, symmetries are used at each level to obtain equivalent 1D models. Section 4 is devoted to the presentation of the discretization method. The instantaneous power conservative structure is first discretized in order to hold the instantaneous power conservation in the discretized model. The discretization of the constitutive relations (derived from the thermodynamics of the mass transfer phenomena) is then conducted in order to guarantee the conservation of energetic properties (conservation or dissipation). Finally Section 5 presents some results of simulation issued from this discretization scheme.

2

Principle of the structured modelling of infinite dimensional thermodynamic systems: a port based approach

In this section we present the port based modelling approach for thermodynamic systems. The structure pointed out in this section is the structure which will be preserved during the reduction procedure proposed in Section 4. We also discuss how the state variables have to be chosen such that geometric properties of the model are emphasized. Indeed, in the Port Based Modelling approach the number of variables is increased in view of pointing out the fundamental geometric power conserving structure of the distributed parameter system, also called Dirac structure [7]. The originality of our approach is to use a geometric power conserving structure for the modelling of dissipative systems. Precisely, in this paper, the geometric structure is used for the modelling of a process where take place convection and diffusion phenomena. For the sake of simplicity, we restrict the presentation to the conservation of matter over a element of volume Ω through which matter is flowing and some reaction or some distributed production of matter (see the example) is occurring. Furthermore we suppose isothermal and isobaric operating conditions. From a mathematical point of view we use the framework of differential geometry [11, 23] to obtain generic coordinate free model that can be used in the case of 2-D or 3-D convection-diffusion processes. More precisely we use differential forms over the spatial domain to define the manipulated variables and the external derivative of these forms to express the driving forces. Differential forms will clarify the relationship between the various quantities that will appear in the paper without be confused by any coordinate system. The main definitions and properties of these mathematical concepts can be found in Appendix A-1. As a consequence, we consider the material balance equation for the species i: Z Z Z ∂qi =− divNi + Nri (1) Ω Ω Ω ∂t 2

where the concentration qi is the state variable, Ni is the flux density and the distributed source flux Nri represents the reaction or adsorption term. qi , N and Nri are respectively a 3-form that can be evaluated on Ω, i.e. qi ∈ Λ3 (Ω), a 2-form that can be evaluated the boundary of the spatial domain ∂Ω, i.e. N ∈ Λ2 (Ω), a 3-form that can be evaluated on Ω, i.e. Nri ∈ Λ3 (Ω). Since the proposed discretization method aims at preserving the structure associated with both material and energy balances, the choice of the manipulated variables concerns moles and energy. First let consider a system in which the thermodynamic variables are spatially uniform and consider the internal energy of a mixture with nc components [36]. It is a function of the extensive variables S, V, ni respectively the entropy, the volume and the particle number. The corresponding intensive variables are the temperature T , the pressure P and the chemical potentials µi . These variables are also function of S, V, ni . They are related by the so-called Euler equation for thermodynamics: U = TS − PV +

nc X

µi ni

U being homogeneous of degree one respect to extensive variables, the variational of U can be written: nc X dU = T dS − P dV + µi dni For simplicity, let consider the Gibbs free energy G = U −T S +P V . G is defined as a function of P, T, ni . The system being isothermal and isobaric, T and P are considered as fixed parameters and the variational of the free Gibbs energy can be written: dG =

nc X

µi dni

(2)

In the distributed parameters case one has to consider the free Gibbs energy density (which is a 3-form), denoted g. From the aforementioned reasons, in the isobaric and isothermal case, g T is a function of the vector of concentrations q = [q1 · · · qi · · · qnc ] i.e. g = g(q) only. Let us R now consider the time derivative of G = Ω g form1 (2): d dt

Z

Z g=

Ω

Ω

∂g T ∂q ∧ = ∂q ∂t

Z µT ∧ q˙ = Ω

Z X Nc Ω

µi ∧ q˙i as

i

∂g =µ ∂q

(3)

Replacing q˙ by its expression in (1) and applying again the Stokes theorem, the time derivative of the free Gibbs energy can be written2 : Z Z d g= µT ∧ (−dN + Nr ) dt Ω Ω Z Z Z (4) =− µT ∧ N + d(µT ) ∧ N + µT ∧ Nr ∂Ω

Ω

Ω

The vector of chemical potentials µ(q) is a vector of 0-forms (functions). If only isothermal and isobaric material balances are considered with null source term i.e. Nr = 0, this leads to the following equation: Z Z Z T T µ ∧N = d(µ ) ∧ N + µT ∧ dN (5) ∂Ω

Ω

Ω

µ and dN being power conjugate variables as well as dµ and N . Equation (5) provides a relation between the internal power variation (right hand side term) and the power flow at the boundary (left hand side term). The relation (5) imposes the choice of the manipulated variables (effort and flow variables) and makes appear the geometric structure representative of the infinite dimensional mass balance, that is to say the Dirac structure that is geometrically defined in [24, 37].As the Dirac structure is the central interconnection structure of our model we are now 1 the free Gibbs energy g = µT ∧ N is considered taking all the species into account and (δ g) and q˙ are q vectors of differential forms 2 We use here the antiderivative property of d say d([ω , ω ]) = [dω , ω ] + (−1)l [ω , dω ] for integration by k l k l k l parts.

3

going to define it precisely for a simple case. Indeed for simplicity reasons we consider the 1-D case using some symmetries (that will be specified in Section 3.2). That means that Ω is assimilated to Z = [0, L] and 3-forms to 1-forms. From the aforementioned reasons the chosen variables for the structured modelling are the following: µ ¶ µ ¶ f1 dN • The 1-forms = are made up with the flow density (extensive variable) e2 dµ and the chemical potential gradient (intensive variable). µ ¶ µ ¶ e1 µ • The 0-forms = are made up with the chemical potential (intensive varif2 N able) and the flux density (extensive variable). Consequently we obtain:

µ

f1 e2

¶

µ =

0 d d 0

¶µ

e1 f2

¶ (6)

Let us define by (f∂ , e∂ ) the boundary power conjugated port variables as a linear combination of the boundary conditions µ∂ and N∂ ∈ ∧0 (∂Ω) such that Stoke’s theorem (5) is satisfied [13], says: Z Z f1 e1 + e2 f2 = e∂ f∂ Ω

Ω

µµ

That means that instantaneous power is null. The quadruple that :

µ

f1 e2

¶

µ =

0 d

d 0

¶µ

f2 e1

¶

µ and

f∂ e∂

f2 e1

µ

¶ , e∂ , 

¶

=

−e1 ∂

f1 e2  

¶

¶ , f∂

such

(7)

f2 ∂

defines a Dirac structure D [12]. This geometric power preserving structure links the pairs of port variables and relates the distributed aspect of the modelling. It has to be completed with two closure equations between the effort and flow variables related with accumulation and dissipation terms. The resulting structured model is depicted on Figure 1.

C

µ

µ

Nr

e1 = µ

0 dq dt

D

1 f 2= N

f1 =dΝ N6

-dµ

e2=d µ

R

N

µ6

Figure 1: Structured representation of the transport phenomena.

The left hand part of the figure named C (according to the electrical circuit notation) is related to the accumulation and allows to compute q from f1 and Nr by time integration and e1 = µ from q using a thermodynamic constitutive equation (first closure equation). The right hand side of the figure named R (according to the resistive element) represents the dissipation term (diffusion or dispersion) and allows to compute f2 from e2 using the second phenomenological law (second closure equation). The D (Dirac) element represents the Dirac structure (interconnection structure). The junctions between C and D and between D and R express the continuity of effort (denoted 0 in figure) or flux variables (denoted 1 in figure) Let us now present the physical system on which will be applied the structured modelling and the associated discretization.

3 3.1

Application to an adsorption process General description

The adsorption process is based on the ability of a solid to preferentially adsorbs the molecules of some species present in a fluid phase and consequently to separate theses species. The 4

central part of the plant associated with the separation process by adsorption is constituted by columns (cylinders) filled with adsorbent pellets (see [15] and [34] for the general presentation of adsorption). Let us consider that the adsorbent is a zeolite and a binary gas mixture where one constituent, named A, is adsorbed faster than the constituent B. The reasons of this selectivity are due on one side to the sorption thermodynamic equilibrium, and on the other side to the kinetic of the mass transfer. Afterward, the gas which is not adsorbed can go through the column ”freely” whenever the other gas is captured and goes to saturate the adsorbent. When all the adsorbent is saturated with the adsorbed gas, one has to stop this operation since the two gases are no more separated. For a detailed description of such processes the reader is referred to [34, 35]. In the sequel we shall consider an isothermal adsorption column and provide a model for the mass transfer phenomena in this system. The mass transfer phenomena description for a bidisperse pellet may be decomposed at three different scales : namely the column scale, the pellet scale and the crystal scale(Cf. Figure 2). Extra granular phase L¼25 cm R int ¼1 cm

rp

z Bidisperse pellet

Macropore scale R p ¼1,24 mm

r Micropore scale Rmic ¼1 µm

Microporous crystal

Figure 2: The three different scales of a bisdisperse pellet. This is the classical approach in Chemical Engineering, justified by the physical sizes of the objects (for instance according to [8], the radius of a crystal is of an order of magnitude of 1 µm and the radius of a pellet approximately 0.8 mm ). For the sake of simplicity we shall restrict our attention to the smaller scale consisting of the microporous medium. The generalization to other scales can be done using the same approach. Yet the final simulation is done on the complete model of the column. Furthermore, we consider without loss of generality that only one species diffuses on the crystal phase. As a consequence the subscript i will be omitted.

3.2

Structured modelling of the microporous scale

We showed in Section 2 that the intrinsic structure of the model and the power exchanges can be pointed out choosing the appropriate set of variables to describe the system. Considering the microporous scale phenomena defined on the spatial domain V mic the chosen variables are the following: • N mic stands for the molar flux of the considered species within the micropore. N mic is a 2-form on the spatial domain V mic ∈ R3 . • µmic stands for the chemical potential and is a 0-form on V mic . • q stands for the concentration within the micropore and is a 3-form on the spatial domain V mic . From homogeneity properties, mass balance is expressed in spherical coordinates (r, θ, φ) and spherical symmetry is used to reduce R3 to R. The spatial domain is then noted R = [0, Rmic ]. As a consequence : • The molar flux N mic becomes a 0-form ∈ ∧0 (R) defined as f2mic = 4πr2 N mic . • The chemical potential (0-form) µmic remains 0-form ∈ ∧0 (R) denoted emic on R. 1 • The concentration initially considered as 3-form, becomes 1-form ∈ ∧1 (R) denoted q L = 4πr2 q. For completing the structured modelling of the crystal as on Figure 1 we also need to define the power conjugate variables: 5

• f1mic (power conjugated to emic 1 ) which states for the divergence of flux and which is a 1-form. • emic (power conjugated to f2mic ) which states for the gradient of the chemical potential 2 and which is defined as a 1-form. Considering the chosen coordinates, the exterior derivative of the 0-form f2mic becomes: ∂f2mic dr ∂r

df2mic =

mic mic One can easily check that the products hf1mic , emic , e2 iL2 [0,Rmic ] are 1 iL2 [0,Rmic ] and hf2 homogeneous and have the dimension of power. We are now going to detail the structure of the microporous scale depicted on the Figure 3.

emic 1b

C

emic 1

emic 1 0 -f

mic 1

f

mic 1

f mic 2b

D

-emic 2 1

f

emic 1a (2) Accumulation

emic 2 mic 2

R

f mic 2

f mic 2a

(1) Interconnection

(3) Dissipation

Figure 3: Structure of the microporous scale model. (1) The interconnection structure - From the previous variables definition one can write : ∂f2mic ∂emic 1 (8) dr, emic = dr 2 ∂r ∂r which is the set of the basic relations constitutive to the Dirac structure defined in Section 2. These equations have to be completed by defining the boundary port variables: f1mic =

mic f∂ = −f2|∂

e∂ = emic 1|∂

(9)

in such a way that Stoke’s theorem can be applied on the crystal volume. These two sets of equations define the interconnection structure and the associated structured model. The mass balance equation and the diffusion equation will form the remaining two closure equations that will complete this model. (2) Closure equation associated with the balance equation and thermodynamic equilibrium - The conserved variable is the linear concentration, the 1-form q L ∈ Λ1 (Z) which obeys the conservation equation: ∂q L ∂f mic = − 2 dr = −f1mic (10) ∂t ∂r The thermodynamic properties give the relations between q L and emic 1 . The Langmuir model [20] is chosen to represent the thermodynamic equilibrium in the adsorbed scale. It links q L to emic by: 1 ¶ µ 1 ∗q L mic 0 (11) e1 = µ (T, P0 ) + RT ln P0 (k ∗ qsL − k ∗ q L ) This equation expresses the chemical potential of the adsorbed species at some temperature T and pressure P0 . k is the Langmuir constant and qsL is the saturation concentration. Remark 3.1 The Hodge star operator ∗ is a linear operator mapping p forms on V to (n − p) forms i.e. : ∗ : Ωp (V ) → Ωn−p (V ) In cartesian coordinates, consider the functions g(z) and the 1-form f (z) = g(z)dz then ∗f (z) = g(z). qsL and q L being 1-forms, ∗qsL and ∗q L are 0-forms. 6

(3) Closure equation associated with the diffusion - We use Knudsen law [19] to present the diffusion in the adsorbed phase (microporous scale). That is to say we only consider the friction exerted by the solid on the adsorbed species. So the constitutive relation representing the diffusion is: f2mic = −

Dmic ∗ q L ∗ emic 2 RT

(12)

where Dmic is the diffusion constant.

4

Model reduction based on geometrical properties

4.1

General concept

As previously mentioned, the port based modelling of the adsorption column is based on basic elements having well defined energetic behavior as already depicted in Figure 3. The distributed aspect of this port based model is essentially supported by the Dirac structure which links power exchanges within the spatial domain and through the boundaries. The proposed discretization method consists in splitting the initial structured infinite dimensional model into N finite dimensional sub-models (finite elements) with the same energetic behavior (cf Figure 4). Furthermore, the support functions used for the interpolation of both effort and flow variables are different to have enough degrees of freedom to guarantee the conservation of the structural properties.

Rmic Power flow exchanged with the following volume

eb1

b b a

C

0

-f

a

0

eab 1

eab 1 ab 1

f ab 1 ea1

f b2

D

eab 2 ab

-eab 2 1

f ab 2

R

ab

f ab 2

f a2

Power flow exchanged with the previous volume

Figure 4: Principle of the spatial discretization

The interconnection between these sub-models is done using the power conjugate boundary port variables. They correspond to the energy flowing across the boundary of one submodel to the boundary of the next sub-model. To each submodel is associated the same generic structure (and consequently parts of the global mass and energy balances) as the global structure, the difference lying in the fact that submodels are finite dimensional i.e. the Dirac structure Dab on Figure 4 is finite dimensional and the reduced effort and flow variables (eab , f ab ) are no more spatially distributed. For simplicity reasons superscript mic will be omitted in the remaining of the section. Let us now explain the reduction scheme.

4.2

Approximation of flows and efforts

The approximation method is based on the separation of variables method. The chosen approximation bases are different according to the degree of the considered differential forms. The 1-form variables f1 and e2 are approximated on Rab = [a, b] ⊂ R = [0, Rmic ] by: f1 (t, z) = f1ab (t)wfab1 (r)

and

7

ab e2 (t, z) = eab 2 (t)we2 (r)

(13)

where the support 1-forms wfab1 and weab2 are chosen such that: Z Z wfab1 = weab2 = 1 Rab

(14)

Rab

R R which implies that f1ab (t) = Rab f1 (t, z) and eab 2 (t) = Rab e2 (t, z). The 0-form variables e1 and f2 are approximated on Rab by: ½ e1 (t, r) = e1 (a)wea1 (r) + e1 (b)web1 (r), (15) f2 (t, r) = f2 (a)wfa2 (r) + f2 (b)wfb2 (r) where the support 0-forms wea1 ,web1 ,wfa2 and wfb2 are chosen such that: ½ a we1 (a) = 1, wea1 (b) = 0, web1 (a) = 0, web1 (b) = 1 wfa2 (a) = 1, wfa2 (b) = 0, wfb2 (a) = 0, wfb2 (b) = 1,

(16)

This choice is done such that e1 (t, a) = e1 (a) and e1 (t, b) = e1 (b).

4.3

Discretization of the interconnection structure

The interconnection structure of each submodel is finite dimensional and concerns the reduced ab ab variables f1ab , eab 2 and their power conjugate variables denoted f2 , e1 . It is defined as follows.  ab   ab  e1 f1  f2ab   eab     2  Definition 4.1 Let us consider the variables   ea∂  ,  f∂a  lying in the discretized eb∂ f∂b version of the Dirac structure (7), then there exist two matrices E and F satisfying (see [12]):  ab   ab  f1 e1   eab  f2ab    2  (17) E  ea∂  + F  f∂a  = 0. b b f∂ e∂ with [E | F ]

full rank

and

EF T + F E T = 0

We require that the approximation variables satisfy the relation induced by the interconnection structure (7). Let us consider the relations induced by the interconnection structure (7) and the approximations (13) and (15):  h i   f 1 = f1ab (t)wfab1 (r) = df 2 = d f2 (a)wfa2 (r) + f2 (b)wfb2 (r) (18)  £ ¤  ab a b e2 = eab (t)w (r) = de = d e (a)w (r) + e (b)w (r) 1 1 1 e2 e1 e1 2 Integrating (18) leads to the following relations between the approximated variables f1ab , eab 2 and e1 (a), e1 (b), f2 (a), f2 (b) :  ab  f1 = f2 (b) − f2 (a) = f∂b (t) − f∂a (t) (19)  ab e2 = e1 (b) − e1 (a) = −eb∂ (t) − ea∂ (t) Equation (19) gives the relations between f1ab , eab 2 (i.e. integrals of the 1-forms f1 , e2 on the spatial domain) and the port variables. Equations (18) and (19) imply that the support functions have to satisfy: dwfa2 (r) = −wfab1 (r), dwfb2 (r) = wfab1 (r) dwea1 (r) = −weab2 (r), dweb1 (r) = weab2 (r)

(20)

We are now going to express the power conjugate variables f2ab , eab 1 (average values of the zero forms on the spatial domain) such that the power on Rab can be written relatively to the finite dimensional variables of the Dirac structure as: £ b b ¤ ab ab ab a a P ab = eab (21) 1 f 1 + e2 f 2 + e∂ f ∂ − e∂ f ∂ 8

Let us write this net power using the previous approximation formulae: Z Z P ab = e1 (r)f1 (r) + e2 (r)f2 (r) + [eb∂ f∂b − ea∂ f∂a ] ab ab Z = {e1 (a)wea1 (r) + e1 (b)web1 (r)}f1ab (t)wfab1 (r) ab Z ab b b a a + {f2 (a)wfa2 (r) + f2 (b)wfb2 (r)}eab 2 (t)we2 (r) + [e∂ f∂ − e∂ f∂ ] ab Z Z = f1ab (t){e1 (a) wea1 (r)wfab1 (r) + e1 (b) web1 (r)wfab1 (r)} ab ab Z Z ab a ab + e2 (t){f2 (a) wf2 (r)we2 (r) + f2 (b) wfb2 (r)weab2 (r)} + [eb∂ f∂b − ea∂ f∂a ] ab

(22)

ab

R

Let us note αab = ab wea1 (r)wfab1 (r). From the compatibility conditions between the supR R port functions (20) one can check that ab wfa2 (r)weab2 (r) = ab web1 (r)wfab1 (r) = 1 − αab and R wb (r)weab2 (r) = αab . Consequently : ab f2 eab 1 ab f2

= [αab e1 (a) + (1 − αab )e1 (b)] = αab ea∂ (t) − (1 − αab )eb∂ (t) = (1 − αab )f2 (a) + αab f2 (b) = (1 − αab )f∂a (t) + αab f∂b (t)

(23) (24)

The relations (19), (23) and (24) linking the ”average” (on the considered spatial domain) ab ab ab values eab 1 , f1 , e2 , f2 to the boundary variables can be summarized as follows: 

1  0   0 0 |

0 −1 0 0

−αab 0 0 1 {z

  ab e1 1 − αab   f2ab 0  a   e∂ 0 1 eb∂ }





0   0 +   1 0 |

0 0 0 1

0 1 − αab 1 0 {z

E

  ab f1 0  eab αab   2 −1   f∂a 0 f∂b }

  =0 

(25)

F

One can easily check that the previous relation relates with the definition of finite dimensional Dirac structure (cf Equation (17) of Definition 4.1) . 4.3.1

Discretization of the diffusion equation

The purpose of this section is to compute f2ab from eab 2 such that energy behavior of the element Rb ab = e ∧ f . We showed, (cf. Equation 12) that from Knudsen law, f R is preserved, i.e. eab 2 2 2 a 2 the diffusion can be written as3 : f2 = −

D ∗ qL D ∗ qL ∗ e2 = −R ∗ e2 , R = RT RT

(26)

Let us now consider the instantaneous power of this element. It can be written in approximated form : Z ¡ ab ¢¡ ¢ ab ab ab 2 PR = e2 (t)weab2 (r) −Reab (27) 2 (t) ∗ we2 (r) = −Kab R(e2 (t)) ab

with

Z Kab = ab

∗weab2 (r)weab2 (r)

After identification the average value f2ab can be written: ab

f2ab = 3e

2

∂P R = −2Kab Reab 2 (t) ∂eab 2

being a one form on the 1D domain R it can be written e2 =

ab ab In the approximated form ∗e2 = eab 2 (t) ∗ we2 (z) = e2 (t)

9

ab ∂we (z) 2

∂r

∂e2 dr ∂r

(28) and ∗e2 =

∂e2 . ∂r

Similarly ∗q L =

∂q L . ∂r

4.3.2

Discretization of the accumulation

The principle is the same here as for the dispersion element. The purpose is to compute eab 1 from q L . The energy on [a, b] of the accumulation element can be written: ¶ Z t Z t µZ Hab = hq˙L , e1 idt = q˙L e1 dr dt (29) 0

0

ab

Let us note that the linear concentration q L lies in the same space as the flow variable q˙L and consequently can be written: q L (t, r) = qab (t)wfab1 (r) (30) The linear saturation concentration qsL (t, r) in the solid phase can be written using the same support function: ab qL (31) s (t, r) = qab,s (t)wf1 (r) From the previous interpolation formulae : µ ¶¸ Z t Z · 1 ∗q L 0 H ab = q˙ab µ (T, P0 ) + RT ln ) wfab1 (r)dt L L P (k ∗ q − k ∗ q ) 0 0 ab s " à ! #Z Z t a qab ∗ wfab1 1 0 = q˙ab µ (T, P0 ) + RT ln ) wfab1 (r)dt P0 (kqab,s ∗ wfab1 − kqab ∗ wfab1 ) 0 b · µ ¶¸ Z t qab 1 0 = q˙ab µ (T, P0 ) + RT ln ) dt P0 (kqab,s − kqab ) 0 {z } | Z = 0

t

(32)

eab 1

q˙ab eab 1 dt

L Consequently eab 1 is deduced from q by: ¶¸ · µ ∗q L 1 ab 0 e1 = µ (T, P0 ) + RT ln P0 (kqab,s − kqab )

With respect to time derivation we obtain : dH ab dH ab = q˙ab eab = eab 1 ⇒ 1 dt dqab

5

(33)

Reduced model and simulation

Relations between the internal average variables and the port variables have been established in the previous section. The balance and closure equations have been established on these reduced variables such that the total mass balance and the energetic behavior of each element are conserved. It is now possible from the boundary conditions to compute the dynamic evolution of each internal variables.

5.1

Simulation procedure

The initial infinite dimensional structured model of the crystal has been discretized taking into account the preservation of its energetic structure. Suppose that the spatial domain has been divided into N sub-domains. On each of the sub-domains has been developed a finite dimensional model with the same structure as the initial model (Figure 5). Concerning the crystal the boundary conditions are the following: • The source flux at the center of the crystal is equal to zero, ı.e. f2 (0) = 0. • The chemical potential at the surface of the crystal is imposed by the macroporous level, ı.e. e1 (Rmic ) = µin . The purpose of the following design procedure is to compute the composition of the adsorbed species along the crystal i.e. the q ab on each segment [a, b]. 10

Surface of the Crystal

Rmic

en1

Boundary Condition

C

en-1n 1

en-1n 1 0

f n-1n 1

-f n-1n 1

C a

C

D

ab

f a2

e11

f 12

D

0 f 01 1 e01

01

-eab 2 1

f ab 2

f ab 2

e01 2

-e01 2 1

f 01 2

f 02

R

f n-1n 2

eab 2

ea1

e01 1

-f 01 1

1

f b2

f ab 1

-en-1n 2

f n-1n 2

eb1 0

e01 1

n-1n

f n-1 2,n

eab 1

-f ab 1

D

en-1n 2

en-1 1,n b

eab 1

f n2

R

R

f 01 2

Boundary Condition

Center of the Crystal

0

Figure 5: Principle of the spatial discretization

Procedure 5.1 (Simulation procedure) We first choose the discretization parameters i.e.: • N the number of discretization steps, • the interpolation functions and forms. The simplest choice compatible with the condition (20) and with normalization (16), is e2 – For the 1-forms: ωab =

– For the 0-forms: ωae2 =

1 b−a dr b−r b−a and

ωbe2 =

r−a b−a .

• the value of αab . Let us note that αab = 0.5 corresponds to a centered method. From these parameters we compute Rab = 2Kab R where Z Z 1 1 e2 e2 Kab = ∗ωab ωab = dr = 2 (b − a) (b − a) ab ab We also proceed to the initialization of the average density qab (t = 0). The following steps has then to be proceeded at each time : • From qab , e1ab is computed thanks to the equation (24). • The pairing of port variables are computed from the boundary conditions. In our case the flux at the center of the crystal is set to zero and the effort at the surface is imposed by the macroporous scale. From the effort in Rmic , µin = e1 (b) = −en∂ and e1 ab we compute en−1 using the equation (25): ∂ en−1 = ∂

e1n−1,n + (1 − αab )en∂ αab

We gradually compute all the boundary port efforts i.e. all the ea∂ and eb∂ . • Computation of the flux gradient. From the previous step all the boundary port effort ea∂ , eb∂ are known and from the equation(25) we compute on each Rab : a b eab 2 = e∂ + e∂

11

• Computation of the flux variables using the R element say: f2ab = −2 ∗ Kab Rab eab 2 • Computation of the boundary port flux : f∂n =

f2n−1,n + (1 − αab )f∂n−1 αab

• Temporal integration and computation of e1 ab using the definition of f1ab f1ab = f∂b − f∂a and the relation

5.2

∂eab 1 = −f1ab ∂t

Generalization to the entire multilevel process simulation

In the previous section only the crystal level of the original system has been presented. Yet modelling and simulation are carried out on the entire multilevel process. The other levels (macroporous scale and column scale) are treated using the same approach. Let first recall all the assumptions we made: • For simplicity, we consider a binary mixture constituted of an inert gas and one component that can be adsorbed. This mixture of gases is supposed to be ideal. • The adsorption column is supposed to be at constant temperature and constant pressure. The velocity v of the flowing fluid is also supposed to be constant. • The diffusion onto the surface of the crystal and the diffusion into the macropore volume are represented by using the Knudsen/Maxwell-Stefan formulation [19]. • In the extragranular phase, a dispersion phenomena is taken into account. It is represented with a constant axial dispersion coefficient. • The Langmuir model for the adsorption equilibrium is used. • The column is supposed to be cylindrical. The particles and the crystals are supposed to be spherical. • Cross section of the column is constant and uniform properties of the adsorbent bed throughout the column are considered. • Concentration gradient in the radial direction of the column are supposed to be negligible. • Spherical symmetries are supposed both in the macropore phase and for the adsorbent. As mentioned modelling and discretization method have been applied to the macroporous scale and to the column scale. The main difference with the crystal modelling lies in the fact that more than one component has to be considered and consequently the diffusion model has to be changed. The Knudsen law is replaced by the Maxwell-Stefan equation. Furthermore, one has to take the convection and dispersion phenomena into account for the modelling of transport phenomena into the column scale. It can be considered as a distributed source term added to the diffusion term. It doesn’t modify the model reduction scheme. When the entire model has to be considered, one has to connect the different models corresponding to the three different scales. The assumption for the interconnection of levels is that the crystal repartition into the pellets is uniform (with porosity4 coefficient ²c ) as well as the pellet repartition within the column (with porosity coefficient ²). Furthermore all the boundaries are at local equilibrium with the outside. Considering spherical symmetry, the crystals have the same chemical potential than the pellet at each point. This means that there is an 4 The

porosity is given by the ratio between the occupied volume and the total volume.

12

equality of efforts at the boundary. Continuity of flux at the boundary of the crystal is also considered. As a consequence, the flux coming from all the crystals present in the sub domain can be seen as a distributed source of flux for the pellet. Considering the geometry of the crystals and the homogeneity of crystal repartition within the pellet, the interconnexion of scales can be formulated as a Dirac structure [9]. Due to this structured modelling, the different scales can be interconnected using port variables with physical meaning. The nature of these port variables avoids numerical difficulties traditionally encountered when unconsistent variables are chosen for the interconnection (numerical stiffness). Finally the number of meshes is chosen such that the solution is consistent and the global simulation doesn’t take too much computational time.

6

Simulation and numerical results

a 10 8 6 4 2

0

100

200 Time (s)

300

400

b 10 8 6 4 2 0

0

100

400

0.25

60

Oxygen molar fraction (% )

Concentration: adsorbed scale (mol /m3)

300

d

c 70

50 40 30 20 10 0

200 Time (s)

0

100

200 Time (s)

300

0.2 0.15 0.1 0.05 0

400

1 0.95 0.9 0.85 0.8 0.75

0

100

200 Time (s)

300

Nitrogen molar fraction (% )

0

Concentration : maroporous scale (mol/m3)

Concentration :extragranular scale (mol/m3)

The simulation of this model has been performed for the separation of a mixture of two constituents, O2 and N2 . For this purpose the model has been spatially discretized by the mixed finite element method we proposed in this paper. The model was programmed in Fortran language and the routine DASPG 5 [28] was used for the numerical integration. Each scale of the column is discretized in ten meshes. The simulation is carried out with the physical parameters presented in [1, 2]. simulated experiment is the response of an adsorption column, initially saturated in N2 , to a steam of air. So in the gaseous phase the concentration profiles are initialized with 0.10 mol.m−3 for the first constituent (O2 ) and with 25.34 mol.m−3 for the second (N2 ). In the adsorbed phase the concentration profiles for the adsorbed species (O2 ) are initialized with 16.3 mol.m−3 . With these initial conditions the process is at thermodynamic equilibrium which corresponds to µ = 0 all along the column.

400

Figure 6: Response to a step change of concentration -(a) in the extragranular scale -(b) in the pellet -(c) in the crystal -(d) at the outlet of the column As shown in fig.(6-d), the output of the column is initially saturated with N2 . When time increases, the concentrations of the two components tend towards those corresponding to the composition of the air i.e. 79% of nitrogen and 21% of oxygen. The response to the change of the concentration of O2 in the ten meshes of the extragranular scale is shown in fig.(6-a). In fig.(6-b), the O2 concentrations in the macroporous medium (in each of the ten mesh associated with the pellet) attached to the first, the sixth and the last discretized mesh of the column are 5 DASPG

is a routine based on Petzold-Gear BDF method. It can be found in IMSL Library

13

presented. In Fig.(6-c), we plot the concentration in the microporous medium (in each mesh) attached at boundaries of the pellet (the last mesh) which is himself attached to the first, the sixth and the last discretized mesh of the column. From the numerical point of view, the use of Port Based modelling and the associated discretisation method avoid stiffness associated with interconnection of multiscale subsystems. Indeed, the variables used for the interconnection of subsystems are the boundary port variables which are the appropriated thermodynamic variables. The interconnection is done using equality of chemical potentials and continuity of flux. When another traditional approach is used, the connection is done using concentrations, which are not continuous at the interface. Consequently one has to use linear interpolation corresponding to the well known film model. This interpolation induces numerical stiffness that may considerably slow down the simulations. Moreover, the discretization method tends to preserve the geometric structure of the model and the energy repartition within the system. As a consequence it is natural to think this strategy will give a good tradeoff between low and high frequency modes. A first study of the dynamic behavior of the discretized model can be found in [32]. In [32] is only considered the simple 1-D ∂q extragranular scale with lumped pellets, i.e. equation (1) with N = vq −Da ∂z , Nr = k(Kq −qp ) (Da = 1e−5 , v = 0.1, k = 1.0, K = 0.5, q(0, z) = 0.2, L = 0.8). It is shown that, qualitatively, the structured model reduction gives better results than finite differences method. Indeed, in Figure 7 are represented the concentration responses along the column respect to an input concentration step. Structural method N=10 1.1

1

1

Concentration of component Q: mol/Volume

Concentration of component Q: mol/Volume

Finite difference with N=10 1.1

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2 0.1

0.9 0.8 0.7 0.6 0.5 0.4 0.3 0.2

0

20

40

60

80

0.1

100

Time: sec

0

20

40

60

80

100

Time: sec

Figure 7: Approximated solutions using finite difference (left) and energy based method (right) N=10. In the case of centered finite differences method oscillations are artificially induced by the discretization. It is not the case with the proposed strategy. Finally, the discretization scheme has been shown to converge to the infinite dimensional spectrum in the centered case.

7

Conclusion

In this paper we proposed a unified procedure for the modelling and spatial discretization of distributed parameter irreversible thermodynamic systems. It is based on the geometrical structure associated with mass and energy balances. In order to illustrate this approach, a model of an adsorption column has been derived directly from its thermodynamical description. The modelling methodology presented exhibits some interesting features : • The model is coordinate free. The geometric specifications may be addressed at the end of the physical modelling process as the derived equations and interconnection structure are ”geometry-independent”. • The model is a network model where each element represents a specific phenomenon which may be identified from a thermodynamics point of view. For instance, in our model, constitutive equations of these subsystems are directly related to the Knudesn/Maxwell Stefan law for diffusion, and the Langmuir’s model for potentials. • The instantaneous power conservation and the description of the power transfers within the system and through its boundaries are explicitly represented. These properties of the model have several important consequences : 14

• The derived model requires parameters that have a clear physical meaning. This considerably simplifies the parameters estimation task. • The model is acausal, hence postpones the choice of boundary conditions (for instance depending here on the model of the gaseous phase in the adsorption column) and is thus clearly reusable. • The central geometric Dirac structure is a direct generalization of Poisson structure in Hamiltonian systems. It suggests and allows the use of passivity-based or energy-shaping techniques for control purposes. These considerations strongly encourage the development of a discretization method which preserves both the nature of the interconnection structures and the physical properties of the connected elements. Such a method has been presented in this paper. The main ideas of this method are the following : • The scheme uses two distinct approximation bases for the power-conjugated variables : one for the intensive variables (functions of the spatial coordinates) and one for the extensive variables (differential forms). This is required by the Dirac structure compatibility conditions and seems quite natural. • The approximation scheme is restricted first in order to preserve the interconnection structure geometric characterization. This means that the reduced interconnection structure is forced to be a discrete Dirac structure. • Then connected discretized elements are forced to be of the same thermodynamics nature as their infinite-dimensional counterparts : dissipative or conservative according to the cases. Further works are actually developed to evaluate the quality of the approximation from a frequency and temporal point of view. A first comparison with the traditional finite differences method can be found in [33]. It is shown that the proposed method presents better results from a numerical and complexity point of view. But the key point is that we now possess a reduced model which allows a direct use of the geometric and thermodynamics properties of the PDEs model to develop estimation or control algorithms. Both the model and the discretization method apply for a large class of distributed parameters thermodynamic systems.

Appendix A-1 The basic derivations of vector analysis are the gradient, curl and divergence. Theses operators are a special case of a single operator, the exterior derivative d applied to quantities that are differential forms. A differential form is a quantity that can be integrated. The type of integral called for by a differential form determine its degree. As an example, a function that can be evaluated ate any point of a spatial domain is a 0-form. The form f (x)dx, where f (x) is a function, is integrated over a single interval and so is a 1-form. The form g(x, y)dxdy is integrated over a surface and so is a 2-form. The form h(x, y, z)dxdydz is integrated over a volume and so is a 3-form. The exterior derivative d applied to a i-form will permit to obtain a i + 1-form. Finally there exists the exterior product denoted by a wedge ∧. This wedge is often omitted when the forms appear under integrals. The exterior product of a i-form and a j-form is a i + j-form. Furthermore the exterior product of the k-form α and the l-form β is such that: α ∧ β = (−1)kl β ∧ α

15

Nomenclature Symbol D G g H P0 PR k N Nr q qsL r R Rc T V e f αab µ0 µ Ω d δ ∗ ∧ mic

f i ab

Name Diffusion coefficient Gibbs free energy Gibbs free energy density Energy function Pressure Dissipated power function Langmuir constant Flux density Reaction or adsorption source flux Concentration Linear saturation concentration Radius coordinate Perfect gas constant Crystal radius Temperature Volume Effort variable Flow variable Discretization coefficitent Reference chemical potential Chemical potential Spatial domain External derivative Variational derivative Hodge star operator Wedge product Superscript for micropore variables Used for approximated variables Component i Variable associated with mesh ab

Unit m2 .s−1 J J.m−1 J Pa W P a−1 mol.m−2 .s−1 mol.m−3 .s−1 mol.m−3 mol.m−1 m J.mol−1 .K −1 m K m3 J.mol−1 J.mol−1 -

References [1] A. Baaiu, F. Couenne, Y. Le Gorrec, L. Lef`evre, and M. Tayakout-Fayolle. Bond graph modeling of an adsorption column. In Preprint 5th International Conference on Mathematical Modelling MATHMOD, Vienna, Austria, February 5-7 2006. [2] A. Baaiu, F. Couenne, Y. Le Gorrec, L. Lef`evre, and M. Tayakout-Fayolle. Port based modeling of a multiscale adsorption column. Mathematical and Computer Modelling of Dynamical Systems,Special Issue on ”Modelling of Distributed-parameter Systems for Control Purposes”, To appear. [3] A. Baaiu, F. Couenne, L. Lef`evre, Y. Le Gorrec, and M. Tayakout-Fayolle. Energy based discretization of an adsorption column. In proc. of Adchem, Brazil, April 2006. [4] R.B. Bird, W.E. Stewart, and E.N. Lightfoot. Transport phenomena. John Wiley and Sons, New York, 2002. [5] F. Couenne, C. Jallut, B. Maschke, P.C. Breedveld, and M. Tayakout. Bond graph modeling for chemical reactors. Mathematical and Computer Modelling of Dynamical Systems, 12(23):159–174, June 2006. [6] F. Couenne, C. Jallut, B. Maschke, M. Tayakout, and P. Breedveld. Structured modeling for processes: A thermodynamical network theorycomputers. Chemical Engineering, In Press, Corrected Proof, Available online 22 April 2007. [7] T.J. Courant. Dirac manifolds. Trans. Amer. Math. Soc. 319, pages 631–661, 1990. 16

[8] F.A. Da Silva and A.E. Rodrigues. Propylene/propane separation by vacuum swing adsorption. A.I.Ch.E. Journal, 47(2):341–357, 2001. [9] D. Eberard, L. Lef`evre, and B. Maschke. Multiscale coupling in heterogeneous diffusion processes: a port-based approach. Proc. International Conference PhysCon,Saint Petersburg, Russia, August 24-26, 2005. [10] D. Eberard and B. Maschke. An extension of port hamiltonian systems to irreversible systems. In Proc. Int. Conf. on Nonlinear Systems Theory and Control, Stuttgart,Germany, 2004. [11] H. Flanders. Differential Forms with Applications to the Physical Sciences. Number ISBN 0-486-66169-5. Dover Publications, New York, 1989. [12] G. Golo, V. Talasila, A.Van der Schaft, and B.Maschke. Hamiltonian disctrization of boundary control systems. Automatica, 40:757–771, 2004. [13] Y. Le Gorrec, H. Zwart, and B. Maschke. Dirac structures and boundary control systems associated with skew-symmetric differential operators. SIAM Journal on Control and Optimization, 42,5:1864–1892, 2005. [14] F. Harley. Differential Forms with Applications to the Physical Sciences. Dover Publications, 1989. [15] E. Jolimaitre. Mass transfer and adsorption equilibrium study in MFI zeolites:aplication to the separation and debranched hydrocarbons in silicalite. PhD thesis, University of Lyon, Lyon 1, 1999. [16] D. Karnopp, D. Margolis, and R. Rosenberg. System dynamics. John Wiley and Sons, New York, 2000. [17] R. Krishna. Multicomponent surface diffusion of adsorbed species : a description based on the generalized maxwell-stefan equations. Chemical engineering Science, 45(7):17791791, 1990. [18] R. Krishna. Problems and pitfalls in the use of the fick formulation for intraparticle diffusion. Chemical engineering Science, 48(5):845861, 1993. [19] R. Krishna and J.A.Wesselingh. The maxwell-stefan approach to mass transfer. Chemical engineering Science, 52:861–911, 1997. [20] I. Langmuir. The constitution and fundamental properties of solids and liquids. i. solidadsorption of gases by solids[j]. Journal of the American Chemical Society, 38:2221 2295, 1916. [21] D. Leinekugel le Cocq, M. Tayakout-Fayolle, Y. Le Gorrec, and C. Jallut. A double linear driving force approximation for non-isothermal mass transfer modeling through bi-disperse adsorbents. Chemical Engineering Science, 62(15):4040–4053, August 2007. [22] D. Leinekugel-Le-Cocq. Contribution to simplified dynamic modelling of a Pressure Swing Adsorption process. PhD thesis, University of Lyon1, 2004. [23] J. Lelong-Ferrand. G´eom´etrie diff´erentielle. Number 2-225-53988-X. MASSON et C ie , Editeurs, 120, boul. Saint-Germain, PARIS-V I e , 1963. [24] B. Maschke and A. van der Schaft. Port controlled hamiltonian representation of distributed parameter sytems. In Workshop on modeling and Control of Lagrangian and Hamiltonian Systems, Princeton, USA, 2000. [25] B. Maschke and A. van der Schaft. Canonical interdomain coupling in distributed parameter systems: an extension of the symplectic gyrator. Int. Mechanical Engineering Congress and Exposition , New- York, USA, Nov. 2001. [26] B.M. Maschke, R. Ortega, and A. J. van der Schaft. Energy- based Lyapunov functions for forced Hamiltonian systems with dissipation. IEEE Trans. on Automatic Control, 45(8):1498– 1502, August 2000. 17

[27] B.M. Maschke, R. Ortega, and A.J. van der Schaft. Energy-based lyapunov functions for forced Hamiltonian systems with dissipation. In proc 37rd CDC, CDC’98, Tampa, Dec 1998. [28] Visual Numerics. Imsl, international http://www.vni.com/products/imsl/.

mathematics

and

statistics

library.

[29] Zienkiewicz O.C. and Zhu J.Z. The Finite Element Method : Its Basis and Fundamentals, volume Engineering mathematics. Butterworth-Heinemann, 2005. [30] R. Ortega, A.J. van der Schaft, I. Mareels, and B.Maschke. Putting energy back in control. IEEE Control Systems Magazine, 21(2):18–32, April 2001. [31] R. Ortega, A.J. van der Schaft, B. Maschke, and G. Escobar. Interconnection and damping assignment: passivity-based control of port-controlled Hamiltonian systems. Automatica, 38:585–596, 2002. [32] H. Peng, A. Baaiu, F. Couenne, and Y. Le Gorrec. Energy based discretization of a class of distributed parameters processes comparison with classical approach. In proc. of the 46th IEEE CDC, News Orleans, Dec 12-14, 2007. [33] H. Peng, A. Baaiu, F. Couenne, and Y. Le Gorrec. Energy based discretization of a class of distributed parameters processes - comparison with classical approach. 46th IEEE Conference on Decision and Control, New Orleans, Louisiana USA, December 12-14, 2007. [34] D.M. Ruthven. Principles of adsorption and adsorption processes. John Wiley and Sons, New York, 1984. [35] D.M. Ruthven, S. Farooq, and K.S. Knaebel. Pressure Swing Adsorption. John Wiley and Sons, New York, 1994. [36] S.I. Sandler. Chemical and Engineering Thermodynamics. Third edition, Wiley and Sons, 1999. [37] A. van der Schaft and B.Maschke. Hamiltonian formulation of distributed-parameter systems with boundary energy flow. Journal of Geometry and Physics, 42:166–194, 2002. [38] A. van der Schaft and B.Maschke. Compositional modelling of distributed-parameter systems, volume Chapter 4. Springer, 2004. [39] J. Villadsen. Solution of Differential Equation Models by Polynomial Approximation, volume ISBN-10: 0138222053,ISBN-13: 978-0138222055. Physical and Chemical Engineering Science, Prentice Hall, February 1978.

18