Phase Separation Processes Using Orthogonal

Efficient Reduced Order Dynamic Modeling of Complex Reactive and MultiPhase Separation Processes Using Orthogonal Collocation on Finite Elements Panos Seferlis1,2, Natassa Dalaouti3 and Theodoros Damartzis2 1

Department of Mechanical Engineering, Aristotle University of Thessaloniki, University Camp, P.O. Box 484, 54124 Thessaloniki, Greece 2

Chemical Process Engineering Research Institute, Centre for Research and Technology – Hellas, 6th km Charilaou-Thermi Road, P.O. Box 60361, 57001 Thermi-Thessaloniki, Greece 3

Testing Research & Standards Center, Public Power Corporation S.A., 9 Leontariou Str., 15351, Kantza-Pallinis/Attiki, Greece 1.1 Introduction Distillation and absorption columns are the most common process units in chemical plants. Separation processes especially those that involve chemical reactions and multi-phase systems give rise to complex dynamic models that require significant computational effort for solution. The main objectives for an absorption or distillation dynamic model are the calculation of accurate predictions of the process behavior during dynamic transition and the reliable solution so that online applications can be facilitated. The calculation of accurate predictions is achieved through the mathematical representation of the major physical and chemical phenomena occurring within the separation column. Multi-component diffusion, chemical reactions, phase equilibrium and so forth must be accounted for in a solid and representative way. Nonequilibrium (NEQ) (rate-based) models provide the most comprehensive description of phenomena in industrial columns [1]. Obviously, the level of detail in the description of such phenomena can be determined only based on the ability to perform sensible experiments that would provide sufficient information for an adequate model validation. The derivation of accurate, precise and uncorrelated estimates for the model parameters is essential in order to ensure good predictive power for the developed dynamic process model [2]. In staged separation processes these phenomena based process models are employed in each column stage (i.e. tray). Similarly, packed separation columns can be approximated by an adequate number of discrete equally spaced points within the column resembling a staged analogue [3]. The number of physical stages in a staged column or equivalent stages in a packed column is usually quite large resulting in large scale process models that require substantial computational effort for solution. Such computational effort hinders the broad utilization of the process model in online control and optimization applications (e.g., model predictive controllers, optimization of transition column

trajectories and so forth). Therefore, an effective reduction in the size of the process model is sought that would not jeopardize the model accuracy by preserving the phenomenological elements in the model in a condensed model structure. Several model order reduction methods for nonlinear dynamic separation processes have been developed. Wave propagation theory has been applied to binary and multicomponent distillation [4-6]. The theory identifies patterns in the composition profiles that are approximated by wave equations and give rise to low-order models. The method seems incapable of handling complex phenomena such as chemical reactions and diffusion controlled distillation. Compartmental models [7] divide the column into compartments and select a representative stage within each compartment; namely the sensitive stage. The sensitive stage is assumed to describe the dynamics of the associated compartment. Steady-state balances around envelopes initiating from the condenser and the reboiler for the stripping and rectifying sections, respectively, are required to close the dynamic balances. Aggregated models [8] also divide the column into compartments with the total holdup of the stages included in a compartment assigned to the holdup of a single representative “aggregation” stage of the compartment. The balance equations for the other stages in the compartment are assumed to be in quasi-steady-state. Several researchers applied the method to complex distillation columns for high purity air separations [910]. The numerical performance of distillation aggregated models is thorough analyzed in [11]. The compromise between solution speed and model accuracy through the manipulation of the number of aggregated stages used is carefully judged. Singular perturbation analysis of the tryby-tray column is performed in [12]. A low order nonlinear model is derived for the slow dynamics that govern the overall dynamic behavior. A balanced empirical grammian with Galerkin projection is used for the model reduction of nonlinear systems with application in distillation in [13]. Most of these methods apply the proposed model reduction techniques in columns that utilize an equilibrium process model with Murphree efficiencies accounting for tray nonideality. A comprehensive overview of fundamental and empirical approaches to the dynamic modeling of distillation columns is available in [14]. Orthogonal collocation on finite elements (OCFE) modeling offers a consistent solution within a unified modeling framework for a wide class of complex separation processes [15]. OCFE techniques in staged columns consider the otherwise discrete column domain as continuous, where composition and temperature profiles are continuous functions of position. Gas and liquid phase material and energy balances are satisfied exactly only at pre-selected points within the column (i.e. collocation points). The technique has been initially introduced by Stewart and coworkers [16] but adopted by many other researchers for steady-state and dynamic simulations of complex separation processes [17-19]. OCFE modeling technique can be applied in conventional and reactive, staged or packed absorption and distillation columns [15]. The key advantages are the high quality of approximation with a substantial degree of model order reduction that further enhances computational efficiency.

Section 2 offers a detailed description of the modeling equations for staged and packed, non reactive, reactive, and multi-phase absorption and distillation processes. Section 3 introduces advanced error monitoring procedures during process dynamic transitions for performance enhancement of the NEQ/OCFE process models. Section 4 presents three illustrative case studies of industrial relevance and interest for reactive absorption, reactive distillation and multi-phase reactive distillation units. The case studies not only confirm the merits and strengths but also identify the issues that require special attention in the NEQ/OCFE model formulation. Finally, the epilogue summarizes the conclusions and provides the context for future research in the area of reduced dynamic modeling for separation units. 1.2 NEQ/OCFE Model Formulation The proposed model formulation combines the rigorous non-equilibrium (rate-based), NEQ, modeling equations [1, 20-21] with the model reduction properties of the OCFE technique. A rate-based model variation assumes that all the resistance to mass and heat transfer is concentrated in two thin film regions adjacent to a phase interface as shown in Figure 1.1. Bulk liquid and gas phases are assumed to have uniform properties in all their control volume whereas chemical reactions may occur in both the bulk and the thin film regions. Thermodynamic equilibrium therefore holds only at the phase interface. The interface structure depends on the number of phases present at any given time instance. In multi phase systems, where a second liquid phase is present the maximum number of possible interfaces is three (i.e. gas-liquid I, gasliquid II, liquid I-liquid II) [22]. Such situation implies that each phase is in contact with all other phases as shown in Figure 1.2b. A common situation is that the second liquid phase appears as a dispersed phase within a continuous liquid phase disallowing the contact between the dispersed liquid and the gas phases. In such a case, only two interfaces are possible (i.e. gas-liquid I, liquid I-dispersed liquid II) as shown in Figure 1.2a. G-L interface gas bulk

gas side film

liquid bulk

liquid side film

Figure 1.1: Thin film model representation.

G-L I interface gas bulk

L I-L II interface liquid I bulk

gas side film

liquid I side film G-L interface

liquid II bulk

liquid I side film LI-LII interface

liquid II side film

(a) G-L II interface gas side film liquid II side film liquid II bulk

liquid II side film L I-L II interface

gas bulk liquid I bulk gas side film

liquid I side film

liquid I side film

G-L I interface

(b) Figure 1.2: Thin film model representations for three phase separation systems. (a) Dispersed liquid phase, (b) two continuous liquid phases. The following assumptions are adopted for the development of the non-equilibrium process model: · One-dimensional mass and heat transport across phase interfaces is assumed. · Thermodynamic equilibrium holds only at phase interfaces. · No axial dispersion along the distillation column is considered. · No entrainment of the liquid phases in the vapour phase occurs. · Complete mixing of bulk phases is achieved. · Heat transfer through the thin film layers is negligible. · Pressure drop inside the column is negligible. · For three phase mixtures, no contact between gas and dispersed liquid phase is possible.

Table 1.1 provides a comprehensive list of the process variables involved in the derivation of the models. Table 1.1: Notation. a specific interfacial area (m2/m3) Acol column stage cross section (m2) A reaction activation energy (J/kmol) c molar concentration (kmol/m3) d molar density (kg/m3) Ð Maxwell-Stefan diffusion coefficient (m2/s) G gas molar flow rate (kmol/hr) H stream molar enthalpy (J/kmol) k0 pre-exponential kinetic parameter K K-value in phase equilibrium L liquid molar flow rate (kmol/hr) m component molar holdup (kmol) n number of collocation points in a finite element N molar flux (kmol/(m2s)) NC number of components NE number of finite elements in a column section NP number of time intervals in simulation domain NR number of chemical reactions NT Number of physical stages P pressure (Pa) Q rate of heat loss / heat duty (J/s) r rate of reaction (kmol/(m3s)) R component rate of reaction (kmol/(m3s)) Rg ideal gas constant (8.314 J/(mol K)) s column position coordinate in OCFE formulation t time (hr) T temperature (K) u component internal energy (J) U stream internal energy (J) W Langrange interpolating polynomial x liquid phase mole fraction y gas phase mole fraction Greek letters activity coefficient g

Γ d Δh h ν j Φ Subscripts t Superscripts G Gb Gf GL GLf int L Lb Lf LI LIb LLIf LII LIIb LLIIf LL

thermodynamic factor film thickness (m) stage height (m) film coordinate stoichiometric coefficient phase volumetric holdup fraction (m3 phase/m3 stage) phase split fraction total gas phase gas bulk region gas film in gas-liquid or gas-liquid1 interface gas-liquid interface liquid film in gas-liquid1 interface interface liquid phase liquid bulk phase liquid film in gas-liquid interface first liquid phase (continuous) first liquid phase (continuous) bulk region liquid1 film in liquid-liquid interface second liquid phase (dispersed) second liquid phase (dispersed) bulk region liquid2 film in liquid-liquid interface liquid-liquid interface

1.2.1 Conventional and Reactive Absorption – Distillation Following the OCFΕ formulation as adopted by [18] and later extended to accommodate non-equilibrium reactive distillation [23] and absorption processes [15] the separation column is divided into sections, with each section defined as the part of the column designated by two material or energy streams entering or leaving the column. Each column section is further divided into smaller sub-domains, namely, the finite elements. The number of collocation points within each finite element determines the order of polynomial approximation for the particular column segment. The main feature of the OCFE formulation is that material and energy balances as resulted from the NEQ rate-based equations are satisfied exactly only at the collocation points. The collocation points are chosen as the roots of the discrete Hahn family of orthogonal polynomials [16]. Such a selection ensures that the collocation points coincide with the location

of the actual stages in the limiting case that the order of the polynomial (i.e. number of collocation points) equals the number of actual stages for any given column section. Hence, the OCFE model approaches the complete tray-by-tray model structure as the number of collocation points increases and the two models become identical if the total number of collocation points becomes equal to the number of trays. Lagrange interpolation polynomials are used within each finite element to approximate the liquid and vapor component flow rates, as well as the liquid and vapor stream enthalpies, as follows: n ~ ~ Li (s ) = åW jL (s j )Li (s j ) j =0

n +1 ~ ~ Gi (s ) = åW jG (s j )Gi (s j ) j =1

0 £ s £ NT , with s0 = 0, i = 1,...., NC

(1.1)

1 £ s £ NT + 1, with sn +1 = NT + 1, i = 1,...., NC

(1.2)

n ~ ~ ~ ~ Lt (s )H L (s ) = åW jL (s j )Lt (s j )H L (s j ) j =0

n +1 ~ ~ ~ ~ Gt (s )H G (s ) = åW jG (s j )Gt (s j )H G (s j ) j =1

0 £ s £ NT , with s0 = 0

1£ s £ NT + 1, with sn +1 = NT + 1

(1.3)

(1.4)

where NC is the number of the components, NT is the number of actual stages included within a ~ given finite element, n is the number of collocation points selection in a given element, Li (s ) and ~ ~ ~ Gi (s ) are the component molar flow rates, Lt (s ) and Gt (s ) are the total liquid and gas stream ~ ~ flow rates, and H L (s ) , H G (s ) are the liquid and gas stream molar enthalpies, respectively. The tilde on the symbols denotes approximation variables. Boundary points of the finite element, s0=0 and sn+1=NT+1, are interpolation points for the liquid and vapor phase approximation schemes, respectively. Functions WL and WG represent Lagrange interpolation polynomials of order n+1 given by the expressions: n

W jL (s ) = Õ k =0 k¹ j

s - sk s j - sk

j = 0,...,n

(1.5)

n +1

W jG (s ) = Õ k =1 k¹ j

s - sk s j - sk

j = 1,..., n + 1

(1.6)

Lagrange polynomials W jL (s ) and W jG (s ) are equal to zero at collocation points sk when k ≠ j and to unity when k=j. Within each finite element, a Lagrange interpolation polynomial of different order can be used. The shape and characteristics of the approximated variable profiles (e.g., linear or irregularly shaped profiles, steep fronts) basically determine the required order of polynomial approximation. Stages that are connected to mass or energy streams entering (e.g., feed streams, streams from a pump-around) or leaving (e.g., product draw streams) the column are treated as discrete equilibrium stages so that the effects of discontinuities in the mass and energy flows within the column sections do not influence the continuity and smoothness of the interpolating polynomial schemes within the elements (equations 1.1-1.4). Such a requirement can however be relaxed when the side material and/or energy streams are distributed uniformly along the length of a finite element or column section (e.g., similar mass and/or heat flows to successive stages within an element) [15]. Hence, the effect of the side streams on the column profiles uniformly spreads throughout the entire element domain and does not diminish the quality of the polynomial approximation. The conditions at the element boundaries obey zero-order continuity, which is a plausible assumption for staged units with discrete profiles. However, certain smoothness conditions for the concentration and enthalpy profiles may be employed in the case of packed columns (e.g., first derivative continuity at element boundary). Assuming zero-order continuity, the boundary conditions connecting two consecutive elements yield the following relations: ~ ~ åW (NT )L (s ) = L (s nk

L j ,k

j =0

k

i

j ,k

i

0 ,k +1

)

i = 1,K NC , k = 1,K NE - 1

nk +1 +1 ~ ~ Gi snk +1,k = å W jG,k +1 (1)Gi (s j ,k +1 )

(

)

i = 1,K NC , k = 1,K NE - 1

(1.7)

(1.8)

j =1

~ ~ ~ åW (NT )L (s )H (s ) = L (s nk

L j ,k

j =0

L

k

t

j ,k

j ,k

t

0, k +1

)H~ (s L

0 ,k +1

)

nk +1 +1 ~ ~ ~ ~ Gt snk +1,k H G snk +1,k = å W jG,k +1 (1)Gt (s j ,k +1 )H G (s j ,k +1 )

(

) (

)

j =1

k = 1,K NE - 1

(1.9)

k = 1,K NE - 1

(1.10)

In equations (1.7)-(1.10) the first index for symbol s refers to the corresponding collocation point and the second one to the corresponding finite element whereas nk is the number of collocation points in the kth element. According to equation (1.7), the extrapolated value for the liquid flow rate at the endpoint of the k-th element (i.e. s=NTk) is set equal to the flow rate at the interpolation and boundary point in the following (k+1)-th element (i.e. s0,k+1=0). Similarly, the vapor flow rate at the endpoint of the k-th element (i.e. sn+1,k=NTk+1), which is also an interpolation point for the vapor approximation scheme, is set equal to the extrapolated value for the vapor flow rate at the beginning of the adjacent (k+1)-th element (i.e. s=1). Assuming that uniform hydrodynamic conditions prevail within each finite element, the dynamic mass balances employed at the collocation point sj, are described by the following equations:

( ) = L~ (s

dmLi s j

i

dt

) ( ) ( ( ) ( )

( ) = G~ (s

dmGi s j

i

dt

( ) )

~ - 1 - Li s j + φL s j RiLb s j + NiLb s j αGL Acol Δh

j

j

)

( ) ( ( )

( ) )

( )

~ + 1 - G i s j + φG s j RiGb s j - N iGb s j α GL Acol Δh

i = 1,....,NC ,

i = 1,....NC ,

j = 1,...,n (1.11)

j = 1,..., n

(1.12)

The terms in the left-hand-side of equations (1.11) and (1.12) account for the component molar accumulation (i.e. material hold-up) in the liquid and gas bulk phase, respectively. These are related to the available liquid and vapor volumes corresponding to an equivalent stage by the equations: miL

~ Li s j col sj =φ sj d sj ~ A Δh Lt s j

( )

( ) ( ) (( ))

i =1,...NC , j =1,...,n

(1.13)

miG

(s j )= φ (s j )d ( ) ( ) ( )

i=1,...NC , j =1,...,n

(1.14)

( )

L

L

G

G

~ Gi s j col sj ~ A Δh Gt s j

( )

( )

In the above equations, d L s j , d G s j

( )

φ L s j , φG s j

denote the liquid and gas phase density and

the liquid and gas phase volumetric fraction, respectively, calculated at the ~ ~ ~ ~ conditions prevailing at the collocation points [ Li s j , Gi s j , T s j , P s j ], while Δh stands

( )

( ) ( ) ( )

for the height of the equivalent stage. Total flow rates and mole fractions in the bulk phases are also defined at the collocation points, as follows:

NC

( ) ( )

xiLb

~ Li s j sj =~ Lt s j

( ) ( ) ( )

i =1,...NC , j =1,...n

(1.15)

( ) ( )

yiLb

( ) ( ) ( )

i =1,...NC , j =1,...n

(1.16)

~ ~ åLi s j = Lt s j i =1

NC

~ ~ åGi s j =Gt s j i =1

~ Gi s j sj = ~ Gt s j

The dynamic mass balances in the gas and liquid films, considering the effect of chemical reactions on the mass transfer through the film regions, as well as the boundary conditions connecting the gas and liquid bulk phases with the respective films are also satisfied exactly only at the collocation points. The mass balances in the films can be written in the following form: ¶ciGf (s j ) ¶N iGf (s j ) + - RiGf (s j ) = 0 ¶t ¶η Gf

i = 1,...NC ,

j = 1,...n,

0 < η Gf £ δ Gf (s j )

(1.17)

¶ciLf (s j ) ¶N iLf (s j ) + - RiLf (s j ) = 0 Lf ¶t ¶η

i = 1,...NC ,

j = 1,...n,

0 < η Lf £ δ Lf (s j )

(1.18)

j = 1,...n

(1.19)

The boundary conditions for equations (1.17) and (1.18) are: N iGb (s j ) = N iGf (s j )

N iLf (s j )

η Lf =d Lf ( s j )

yiGb (s j ) = yiGf (s j )

ηGf =0

= N iLb (s j )

xiLf (s j )

ηGf =0

( )

η Lf =d Lf s j

i = 1,...NC ,

= xiLb (s j )

i = 1,...NC ,

j = 1,...n

(1.20)

The diffusion molar flux term Ni(sj) in equations (1.11), (1.12), (1.17) and (1.18) is estimated by the Maxwell-Stefan equations for multicomponent mixtures, at the conditions prevailing at each collocation point. For ideal gas phase and one-dimensional mass transfer normal to the interface, the Maxwell-Stefan equations for the gas and the liquid thin film sides take the following form: Gas phase

¶yiGf (s j ) ¶η

Gf

(y (s )N (s ) - y (s )N (s )) = -å (P~(s )/ R T~(s ))Ð NC

k =1 k ¹i

Gf k

j

Gf i

j

Gf i

j

g

j

Gf k

j

G ik

j

i = 1,...NC,

j = 1....n,

0 < ηGf £ δGf (s j ) (1.21)

Liquid phase NC -1

å G (s ) i ,k

¶xkLf (s j )

j

¶η

k =1

Lf

(x (s )N (s ) - x (s )N (s )) = -å c (s )Ð NC

Lf k

k =1 k ¹i

j

Lf i

Lf i

j

t

j

j

Lf k

j

L ik

i = 1,...NC - 1,

j = 1,...n, 0 < η Lf £ δ Lf (s j )

(1.22) where the gamma thermodynamic factor is defined as follows: ¶ ln g i (s j ) Gi k (s j ) = δ ik + xi (s j ) ¶xk (s j ) ~ ( ) ~ ( ) ( ) T s j , P s j , xk s j ,k ¹i =1,... NC -1

(1.23)

Terms Ri(sj) in Eq (1.11), (1.12), (1.17) and (1.18) denote the total component reaction rate of the ith component in the gas and liquid bulk and film regions and are estimated at the conditions prevailing at each collocation point by the following equations: Ri (s j ) = ån i ,r rr (s j ) NR

i = 1,...NC ,

j = 1,....n

(1.24)

r

where rr ( s j ) denotes the rates of the reactions taking place. OCFE formulation can be tailored to allow reactive and non-reactive sections in the column. At the gas-liquid interface of each collocation point thermodynamic equilibrium is assumed, described by the following equation: yiint (s j ) = K i (s j ) xiint (s j )

i = 1,...NC ,

j = 1,...n

(1.25)

The boundary equations at the interface are: N iGf (s j )

( )

ηGf =d Gf s j

yiGf (s j )

= N iint (s j ) = N iLf (s j )

( )

ηGf =d Gf s j

= yiint (s j )

η Lf =0

i = 1,...NC ,

xiint (s j ) = xiLf (s j )

η Lf =0

j = 1,...n

i = 1,...NC ,

(1.26)

j = 1,...n

(1.27)

Neglecting the heat transfer effects along the film regions the overall dynamic energy balance at each collocation point becomes:

dU(s j ) ~ ~ ~ ~ ~ ~ ~ ~ = Lt (s j -1)H L (s j -1) + Gt (s j +1)H G (s j +1)- Lt (s j )H L (s j ) - Gt (s j )H G (s j ) + Q(s j ) j = 1,...n dt

where U (s j ) = å {miL (s j )u iL (s j ) + miG (s j )uiG (s j )} NC

(1.28)

j = 1,....n

i =1

Term Q(sj) is the net heat transferred from the surroundings. Contrary to the basic principle of the OCFE formulation that isolates stages with a mass or energy streams leaving or entering the column and treats them as discrete stages, the systematic and uniform heat exchange that may take place in a given column section can be also formulated as in equation (1.28). Deviations from the basic rule are only verified when the heat exchange within a given element is distributed uniformly along the element. For instance, heat losses or heat exchange for a group of neighboring stages can be approximated using balance equation (1.28) within finite elements that exclusively contain the stages where the heat exchange takes place. However, pump-around streams that carry significant material and energy amounts from one point of the column to another should be definitely treated as discrete stages. Empirical correlations are employed for the calculation of the pressure drop in the column, the liquid hold-up, the specific interfacial area, and the gas and liquid film thickness. These correlations consider column internals and hydraulics and depend on the type of the column plate (e.g., sieve, bubble cap) or column packing. The calculations are performed for the conditions prevailing at the collocation points. The column is also subject to other operating limitation such as flooding constraints. In conclusion, the NEQ/OCFE model formulation preserves the detailed description of the physical and chemical phenomena occurring within an absorption or distillation column as the modeling set of equations for collocation points and discrete stages are equivalent, but uses fewer points for the calculation of the material and energy balances than the full-order model (e.g., trayby-tray model formulation) as shown in Figure 1.3. Therefore, a more compact representation is available without diminishing the predictive power of the process model.

Tray-by-tray (full-order) model

OCFE (reduced-order) model

Model representation per tray or collocation point

collocation points element boundary

two-phase case

three-phase case

Figure 1.3: Full- and reduced-order model representation. 1.2.2 Multi-phase reactive distillation The appearance of a second liquid phase within the distillation column is imposed by the phase thermodynamic stability. The second liquid phase can be assumed either as a dispersed liquid phase within a continuous liquid phase or a second continuous liquid phase. Therefore, the thin film model representation depends on the assumption about the type of the second liquid phase (Figure 1.2). In the current chapter, the continuous-dispersed model for the two liquid phases is considered as the most likely to occur situation. However, the model is readily expandable to accommodate the more complex three phase scenario. A unique feature of multiphase distillation is the fact that the boundary between the three-phase to the two-phase region in the column introduces a discontinuity in the concentration and the temperature profiles. Appearance or disappearance of a second liquid phase results in a change in the model structure as the balance equations for the second liquid phase and the corresponding interfaces should be incorporated or removed from the model, respectively. A second drawback is that the phase boundary is not known and its location may change significantly during a column dynamic transition. NEQ process models have been developed in [22, 24]. The derivation for the NEQ/OCFE process balances for each phase proceeds as follows [25]:

Liquid I phase: dmiLI (s j )

(

)

~ ~ = LiI (s j - 1) - LiI (s j ) + φ LI (s j )RiLIb (s j ) + N iGLb (s j )a GL - N iLLIb (s j )a LL Acol Δh

dt i = 1,....., NC ,

(1.29)

j = 1,......, n

Liquid II (dispersed liquid) phase: dmiLII (s j )

~ ~ = LiII (s j - 1) - LiII (s j ) + (φ LII (s j )RiLIIb (s j ) + N iLIIb (s j )a LL )Acol Δh ,

dt i = 1,....., NC ,

(1.30)

j = 1,......, n

Gas phase: The gas phase balance is derived from equation (1.12). In the dynamic mass balances (1.29) and (1.30), mi denotes the component molar accumulation, which may be computed from the following relations: miLI

~ LiI (s j ) col (s j ) = φ (s j )d (s j ) ~I A Δh i = 1,...., NC , Lt (s j )

miLII

LI

LI

(s j ) = φ (s j )d LII

LII

j = 1,....n

~ LiII (s j ) col (s j ) ~II A Δh i = 1,...., NC , Lt (s j )

j = 1,....n

(1.31)

(1.32)

where φ stands for the phase volumetric holdup fraction and d stands for the phase molar density at each collocation point. The gas phase holdup is calculated from relation (1.14). The total molar flow rates and molar fractions for each phase are computed from the component flow rates as follows: ~ ~ å LiI (s j ) = LtI (s j ) NC

i =1

~ ~ å LiII (s j ) = LtII (s j ) NC

i =1

xiI ,Lb

~ LiI (s j ) (s j ) = ~II i = 1,...., NC , Lt (s j )

(s j ) = L~iII (s j ) Lt (s j )

j = 1,...., n

(1.33)

~II

xiII ,Lb

i = 1,...., NC ,

j = 1,...., n

(1.34)

Under the assumption of negligible heat transfer through the films, the dynamic energy balance attains the form:

dU (s j )

~ ~ ~ ~ ~ ~ = LtI (s j - 1)H LI (s j - 1) + LtII (s j - 1)H LII (s j - 1) + Gt (s j + 1)H G (s j + 1)

dt ~ ~I ~ ~ ~ ~ - Lt (s j )H LI (s j ) - LtII (s j )H LII (s j ) - Gt (s j )H G (s j ) + Q (s j )

(1.35)

j = 1,......, n

In equation (1.35) Q(sj) denotes the net heat rate exchanged between the column and its ~ surroundings, H the stream molar enthalpy and U(sj) the total molar energy accumulation, which is calculated from the relation:

{

}

U (s j ) = å miLI (s j )uiLI (s j ) + miLII (s j )u iLII (s j ) + miG (s j )uiG (s j ) NC

j = 1,....,n

(1.36)

i =1

The dynamic mass balances in the gas and liquid films at the two possible interfaces (G-LI and LI-LII) along with the boundary conditions are as follows: G-L interface – Gas side film ¶ciGf (s j ) ¶N iGf (s j ) + =0 ¶t ¶η Gf

Boundary conditions N iGb (s j ) = N iGf (s j )

η Gf =0

i = 1,...., NC ,

j = 1,..., n , 0 < η Gf £ δ Gf (s j )

yiGb (s j ) = yiGf (s j )

η Gf =0

i = 1,...., NC ,

(1.37)

j = 1,...., n

(1.38)

G-L interface – Liquid side film

¶ciGLf (s j ) ¶N iGLf (s j ) - RiGLf (s j ) = 0 + GLf ¶t ¶η

Boundary conditions

N iLIb (s j ) = N iGLf (s j )|ηGLf =δ GLf (s ) j

i = 1,...., NC ,

j = 1,..., n , 0 < η GLf £ δ GLf (s j )

xiLIb (s j ) = xiGLf (s j )|ηGLf =δ GLf (s ) j

i = 1,...., NC ,

j = 1,....,n

(1.39)

(1.40)

L-L interface – Liquid I side film

¶ciLLIf (s j ) ¶N iLLIf (s j ) + - RiLLIf (s j ) = 0 ¶t ¶η LLIf

Boundary conditions

N iLIb (s j ) = N iLLIIf (s j )|η LLIf =0

i = 1,...., NC ,

xiLIb (s j ) = xiLLIf (s j )|η LLIf =0

j = 1,..., n , 0 < η LLIf £ δ LLIf (s j )

(1.41)

i = 1,...., NC ,

(1.42)

j = 1,....,n

L-L interface – Liquid II side film

¶ciLLIIf (s j ) ¶N iLLIIf (s j ) + - RiLLIIf (s j ) = 0 ¶t ¶η LLIIf

i = 1,...., NC ,

j = 1,..., n , 0 < η LLIIf £ δ LLIIf (s j )

(1.43) Boundary conditions

N iLIIb (s j ) = N iLLIIf (s j )|η LLIIf =δ LLII (s ) j

xiLIIb (s j ) = xiLLIIf (s j )|η LLIIf =δ LLIIf (s ) j

i = 1,...., NC ,

j = 1,....,n (1.44)

Thermodynamic equilibrium is valid at the phase interfaces (G-LI and LI-LII) and is described by the following equations: yiGL (s j ) = K iGL (s j )xiGL (s j ) xiI ,LL (s j ) = K iLL (s j )xiII ,LL (s j )

i = 1,....NC ,

j = 1,...., n

i = 1,....NC ,

(1.45)

j = 1,...., n

(1.46)

In equation (1.46) xiI ,LL and xiII ,LL represent the molar fraction of the first and second liquid phase at the LL interface respectively. The boundary conditions that arise from the equilibrium equations at the interfaces are: N iGf (s j ) yiGf (s j )

( )

= N iGL (s j ) = N iGLf (s j )

( )

= yiGL (s j ), xiGL (s j ) = xiGLf (s j )

η Gf =δ Gf s j

η Gf =δGf s j

N iLLIf (s j ) xiLLIf (s j )

( )

η LLIf =δ LLIf s j

( )

η LLIf =δ LLIf s j

i = 1,...., NC ,

j = 1,...., n

(1.47)

η GLf = 0

= N iLL (s j ) = N iLLIIf (s j )

η GLf =0

i = 1,...., NC ,

i = 1,...., NC ,

j = 1,...., n

j = 1,...., n

(1.48)

(1.49)

η LLIIf = 0

= xiI ,LL (s j ), xiII ,LL (s j ) = xiLLIIf (s j )

η LLIIf =0

i = 1,...., NC ,

j = 1,...., n

(1.50)

The reaction terms Ri used in the above equations represent the total rate of reaction per component i and are computed from equation (1.24). The component molar flux terms Ni at the collocation points are calculated through the multi-component Maxwell-Stefan relations (equations 1.21 and 1.22 for the gas and liquid phase, respectively).

Equations (1.29)-(1.50) and (1.21)-(1.22) form the NEQ/OCFE model for a three-phase element whereas equations (1.1)-(1.29) form the NEQ/OCFE model for a two-phase element. One major feature of three-phase distillation units is that the liquid phase split regions are not known before a solution of the modeling equations is obtained. Besides, during dynamic transition the phase boundary may shift with time. Unless the correct set of modeling equations is used, thus identifying the phase distribution correctly, the simulated results will not represent accurately the column behavior and occasionally fail to reach a feasible solution. The latter is possible if a three-phase element is used to represent a two-phase region. Therefore, a tracking mechanism for the phase boundary is necessary for improved accuracy. According to the methodology of [26] as later implemented in [25] an element breakpoint is adaptively located at the phase transition boundary. The identification of the phase discontinuity point is performed via a liquid-liquid (LL) flash calculation at the element breakpoints. The LL flash calculation determines whether the formation of the second liquid phase is possible and subsequently specifies whether the element breakpoint becomes a phase transition boundary. The element breakpoint must be allowed to move freely during simulation or optimization studies in order to trace variations in phase allocation. The LL flash at the element breakpoints involves the following additional calculations: LL flash equations: Φ(s0 )xiI ,lb (s0 ) + (1 - Φ(s0 ))xiII ,lb (s0 ) = xi (s0 )

å (x (s ) - x (s )) = 0 NC

i =1

I ,lb i

0

II ,lb i

0

γ iI xiI ,lb (s0 ) = γ iII xiII ,lb (s0 )

i = 1,...., NC

i = 1,...., NC

i = 1,...., NC

(1.51)

(1.52)

(1.53)

where γi is the activity coefficient of each component and Φ is the phase split fraction (Φ=1 for single liquid phase and 0