A Non-equilibrium Approach to Model Flash Dynamics with Interface Transport

TL;DR

This paper develops a non-equilibrium thermodynamics-based modeling framework for liquid-vapor systems, capturing interface transport and stability properties through differential-algebraic equations.

Contribution

It introduces a novel DAE-based model for non-equilibrium flash dynamics incorporating interface transport and stability analysis.

Findings

The model is proven to be of index one.

Equilibrium is unstable in non-isobaric regimes.

Isobaric regimes are numerically shown to be stable.

Abstract

This article presents a modeling framework for a class of multiphase chemical systems based on non-equilibrium thermodynamics. Compartmental modeling is used to establish the dynamic properties of liquid-vapor systems operating far from thermodynamic equilibrium. In addition to the bulk-phase molar/energetic dynamics, interface transport processes yield to algebraic constraints in the model description. The irreversible system is thus written as a system of Differential-Algebraic Equations (DAEs). The non-equilibrium liquid-vapor DAE system is proven to be of index one. A local stability analysis for the model shows that the equilibrium state is unstable for non-isobaric operation regimes, whereas numerical evidence shows that isobaric operation regimes are stable. To extend the stability analysis, internal entropy production for the irreversible flash-drum is presented as a Lyapunov function candidate.

Tables

Table 1. Table 1: Numerical values for the stationary state

$y_1^{\star }=0.6615$,	$y_2^{\star }=0.3385$	$T^{\star }={78.09}^{o}C$,	$V_l^{\star }=0.1m^3$
$x_1^{\star }=0.2764$,	$x_2^{\star }=0.7236$	$C_g^{\star }=34.6874\mathrm{mol}/m^3$,	$F_{\alpha ,V}^{\star }=1m^3/s$

Table 2. Table 2: Dimensionless heat capacity 𝚌 v subscript 𝚌 v \mathtt{c}_{\text{v}} for ideal gases

	Gas	Temperature range
${𝚌}_{\text{v}}=3/2$	Monotonic non-interactive atoms	K
${𝚌}_{\text{v}}=5/2$	Diatomic non-interactive molecules	K
${𝚌}_{\text{v}}=7/2$	Diatomic non-interactive molecules	$T>{10}^3$K

Equations

$$S=\text{S}(U,V,N_1,\dots ,N_c)$$

$$\text{S}(\lambda U,\lambda V,\lambda N_1,\dots ,\lambda N_c)=\lambda \text{S}(U,V,N_1,\dots ,N_c),\forall \lambda >0.$$

$$\frac{\partial S}{\partial U}=\frac{1}{T},\frac{\partial S}{\partial V}=\frac{P}{T},\frac{\partial S}{\partial N_j}=\frac{-{\mu }_j}{T},j\in \{1,\dots ,c\}.$$

$$\text{d}S=\frac{1}{T}dU+\frac{P}{T}dV+\frac{{\mu }_1}{T}d{N}_1+\dots +\frac{{\mu }_c}{T}d{N}_c.$$

$$\frac{dS}{dt}=\frac{1}{T}\frac{dU}{dt}+\frac{P}{T}\frac{dV}{dt}+\frac{{\mu }_1}{T}\frac{d{N}_1}{dt}+\dots +\frac{{\mu }_c}{T}\frac{d{N}_c}{dt}.$$

$$\frac{d{N}_{g,j}}{dt}$$

$$\frac{d{N}_{l,j}}{dt}$$

$$\frac{d{U}_g}{dt}$$

$$\frac{d{U}_l}{dt}$$

$$\frac{d{K}_g}{dt}$$

$$\frac{d{K}_l}{dt}$$

$$V_l={\bar{\text{v}}}_1{N}_{l,1}+\cdots +{\bar{\text{v}}}_c{N}_{l,c},$$

$$V_g=V_o-V_l,$$

$$F_{\alpha ,V,\mathrm{out}}$$

$$Q_{\alpha }={\lambda }_{\alpha }(T_{\alpha ,Q}-T_{\alpha }),$$

$$U_{\alpha }=U_{o,\alpha }(N_{\alpha })+{\mathcal{C}}_{\alpha }(N_{\alpha })(T_{\alpha }-T_{\text{o}}),$$

$$P_{\text{g}}(U_g,V_g,N_g)$$

$$P_l=P_i=P_g(U_g,V_g,N_g).$$

$$[T_i,x_{1,i},y_{1,i},\dots ,x_{c,i},y_{c,i},n{]}^{\text{t}}.$$

$$n_{g,j}$$

$$n_{l,j}$$

$$n_{\alpha }=\sum _{j=1}^c{n}_{\alpha ,j}.$$

$$e_g$$

$$e_l$$

$$e_{\alpha ,K}$$

$$v_{\text{i},\alpha }=\frac{m_{\alpha }}{{\rho }_{\alpha }{\text{A}}_i}.$$

$$k_{g,j}{C}_g(y_{j,i}-y_j)-k_{l,j}{C}_l(x_j-x_{j,i})+(y_j-x_j)n=0,$$

$$\sum _{j=1}^c{n}_j\Delta {\bar{h}}_{\mathrm{vap},j}+{\lambda }_{g,i}(T_g-T_i)-{\lambda }_{l,i}(T_l-T_i)=0,$$

$$0$$

$$\dot{\mathbf{z}}$$

$$0$$

$$\mathbf{z}=[N_{g,1}\dots N_{g,c}{N}_{1,l}\dots N_{l,c}{U}_g{U}_l{K}_g{K}_l{]}^{\text{t}},$$

$$\mathbf{w}=[y_{1,i}\dots y_{c,i}{x}_{1,i}\dots x_{c,i}{T}_{i}n{]}^{\text{t}}.$$

$$\displaystyle\sigma=\sum_{\alpha\in\{g,l\}}\bigg{(}F_{\alpha,S,\mathrm{out}}-F_{\alpha,S,\mathrm{in}}-\frac{Q_{\alpha}}{T_{\alpha,Q}}\bigg{)}+\frac{dS}{dt},$$

$$\frac{dS}{dt}=\sum _{\alpha \in \{g,l\}}\frac{d{S}_{\alpha }}{dt}=\sum _{\alpha \in \{g,l\}}\frac{1}{T_{\alpha }}\frac{d{U}_{\alpha }}{dt}+\frac{P_{\alpha }}{T_{\alpha }}\frac{d{V}_{\alpha }}{dt}+\frac{-{\mathbf{\mu }}_{\alpha }^{\text{t}}}{T_{\alpha }}\frac{d{\mathbf{N}}_{\alpha }}{dt},$$

Full text

Keywords: Non-equilibrium thermodynamics, Multiphase systems, Differential-algebraic systems, Flash-drum dynamics, Lyapunov stability, Entropy production

1 Introduction

Despite its historical and contemporary significance, design, control and operation of multiphase processes have been challenging tasks in process systems engineering over the years (Skogestad, 1997; Taylor and Krishna, 2000). In this article, we explore the possibility of modeling a special class of multiphase systems, the open flash-drum, utilizing a non-equilibrium physics-based perspective. In contrast with traditional equilibrium formulations, non-equilibrium models can keep track of irreversible phenomena such as energy degradation and entropy production (de Groot and Mazur, 1984). Moreover, irreversible schemes have led to useful insights for stability analysis and feedback control design in the chemical process systems literature (Alonso and Ydstie, 1996; Ydstie and Alonso, 1997; Favache and Dochain, 2009; García-Sandoval et al., 2015; Ydstie, 2016). The work presented in this article expands the available modeling techniques and stability theory to include open liquid-vapor systems that operate far from thermodynamic equilibrium.

Dissipative systems theory, originally proposed by Willems (1972) as an extension of (linear) passivity-based analysis, has been established as a fundamental tool for analysis and control design of mechanical and electrical systems. As for chemical systems, dissipative systems theory has received an increasing level of attention after Alonso and Ydstie (Alonso and Ydstie, 1996; Ydstie and Alonso, 1997) developed the concept in detail using the first and second laws of thermodynamics. A decade later, in the article by Favache and Dochain (2009), the possibility of characterizing the continuous stirred tank reactor (CSTR) through energetic and entropic formulations is explored with insightful results. Other worth mentioning references include the work by Hoang et al. (2012), where Lyapunov control laws, based on thermodynamic availability, are proposed for a larger class of non equilibrium CSTRs; and the article by García-Sandoval et al. (2015), where dissipative properties, based on internal entropy production, are established for irreversible CSTRs.

Thermodynamic process systems theory has focused mostly on single-phase systems, also known as simple thermodynamic systems. And it is not until recently that multiphase systems theory analysis has gained an increasing interest from the scientific community.

For simple thermodynamic systems, entropy (resp., internal energy) is a strictly concave (resp., convex) function of the extensive variables (Callen, 1985), as depicted Figure 1 (right) below. Such concavity properties had permitted to assess for stability using Lyapunov theory and dissipative analysis (Favache and Dochain, 2009; García-Sandoval et al., 2015; Hoang et al., 2012). For multiphase systems the concavity is not strict (Gromov and Caines, 2012) as the very existence of separate phases follows as a consequence of the loss of concavity in entropy (Callen, 1985), see Figure 1 (left) below. To characterize a system regardless of its convexity properties, we propose to study liquid-vapor processes based on internal entropy production rather than entropy or internal energy potentials.

Historically, systems theory analysis for multiphase processes can be traced back to the pioneering work of Rosenbrock (1963), who demonstrated that a non-ideal binary distillation column operates at a unique asymptotically stable steady state. Three decades later Rouchon and Creff (1993) developed a stability analysis, based on geometric considerations, for a multicomponent flash-drum. Their approach unfortunately does not extend to multistage process units. In a more recent contribution, Ydstie (2016) developed conditions for the existence of a unique stable steady-state for an adiabatic flash-drum operating on an equilibrium manifold. Looking forward to extend the process systems theory to include dynamic non-equilibrium liquid-vapor units, a modeling framework is presented in our previous contribution Romo-Hernandez et al. (2018). The article discussed here serves as an extension to Romo-Hernandez et al. (2018) that includes the construction of the entropy production function for the irreversible flash-drum. As a Lyapunov function candidate, entropy production can serve to characterize the dynamic properties for nonlinear irreversible flash-drums. In addition, details on mathematical procedures excluded from Romo-Hernandez et al. (2018) are included thorough different sections in the work presented below.

In this article, we consider a non-equilibrium flash-drum as the combination of three subsystems: a liquid phase, a gas phase, and an interface. To study such system, the paper is organized as follows. In Section 2, the non-equilibrium liquid-vapor model is built as a system of differential-algebraic equations (DAEs) from mass and energy conservation principles. The DAE model takes into account exchange processes through the interface, viewed as constraints on the dynamics of the system. In Section 3, internal entropy production for the flash-drum is written as the sum of products of flows and driving forces. Being positive definite, internal entropy production is a physics-based Lyapunov function candidate to characterize the stability of the irreversible flash-drum. Numerical results are presented in Section 4. First, the DAE model is rewritten and improved for numerical integration using a bijective change of coordinates. Second, a non-ideal water-methanol liquid gas mixture is considered to present dynamic trajectories and results regarding Lyapunov’s first and second methods for stability. Conclusions and future areas for research are discussed in Section 5.

2 Liquid-Vapor Thermodynamic Systems

In this section, we develop a system of nonlinear differential-algebraic equations (DAEs) of index one that describes a flash-drum as the interconnection of two thermodynamic subsystems. We motivate the need for such description through Gibbs equation as a way to describe entropy variations inside an open liquid-vapor system. The obtained DAE model is an abstract representation of the first law of thermodynamics applied to the irreversible flash-drum.

A thermodynamic system is completely defined once its physical properties are determined. For instance, all the physical properties of a closed chemical system with $c$ chemical components can be recovered once a particular set of coordinates $(U,V,N_{1},\ldots,N_{c})$ is known. Such coordinates represent the internal energy, volume, and mole numbers respectively. This is known as the first postulate of thermodynamics and it was proposed by Callen (1985), on the basis of the Gibbs’ school on modern thermodynamics.

The second postulate of thermodynamics defines entropy as a function that is maximized at thermodynamic equilibrium (Callen, 1985). Moreover, such function relates the entropy with the intensive variables (Wightman, 1979).

Definition 1** (Entropy).**

Let $(U,V,N_{1},\ldots,N_{j})\in\mathbb{R}_{>0}^{c+2}$ represent the internal energy, volume, and mole numbers for a closed thermodynamic system, and let $(T,P,\mu_{1},\ldots,\mu_{c})$ stand as the temperature, pressure and the chemical potential. Then, entropy is a concave function

[TABLE]

As a consequence of (1b), the coordinates $(T,P,\mu_{1},\ldots,\mu_{c})$ are said to be the conjugates to the extensive variables $(U,V,N_{1},\ldots,N_{c})$. $\Box$

The formal structure of thermodynamics, Equation (1), is defined for closed systems at thermodynamic equilibrium (Callen, 1985). Nevertheless, Equation (1) is still considered valid locally when studying systems that are not at thermodynamic equilibrium (Hoang et al., 2012). Non-homogeneous systems, non-stationary processes and open vessels exchanging mass and energy with the environment are some examples of non-equilibrium systems. Computing the differential of (1a) we obtain what is known as the Gibbs equation

[TABLE]

When time variations are considered instead of the differentials in (2) we get

[TABLE]

As the derivatives on the right hand side of (3) can represent exchange rates between a system and its environment, Equation (3) is frequently considered as an extension of Gibbs equation (2) for open thermodynamic systems. In the following section, we discuss how to take advantage of conservation laws to describe the exchange rates for $(U,V,N_{1},\ldots,N_{c})$ that appear on the right hand side of (3), particularly for the open flash-drum.

2.1 Conservation Principles

From here on, we consider the system of interest as a flash-drum with one liquid phase on the bottom of a rigid vessel and a vapor phase above as depicted in Figure 2. Inside the system, there are $c$ chemical components distributed between phases. To avoid unnecessary repetition, we set sub-index $\alpha\in\{l,g\}$ as the phase sub-index.

Each phase inside the flash-drum has certain amount of internal energy $U_{\alpha}$, and a determined number of moles $N_{\alpha,j}$ for each component $j\in\{1,\ldots,c\}$. As the system is open, moles flow in and out at convective rates $F_{\alpha,N_{j},\mathrm{in}}$ (mol/sec) and $F_{\alpha,N_{j},\mathrm{out}}$ (mol/sec) respectively. Thermal energy is exchanged between phase $\alpha$ and the environment at a rate $Q_{\alpha}$ (J/sec) through a heat exchanger. Additionally, moles and energy can be exchanged between phases. We set the interface energy exchange rate as $e$ (J/sec) and the interface molar exchange rate as $n$ (mol/sec). These represent the rates at which total energy and total mole numbers are exchanged from the liquid to the gas phase. As a final consideration, we set $K_{\alpha}$ to be the kinetic energy of the flow moving through phase $\alpha$.

For modeling purposes, the following assumptions are considered:

1. A1.

Each phase is perfectly mixed. 2. A2.

Potential energy is constant all over the process. 3. A3.

Flow compressibility and viscous losses are negligible. 4. A4.

Kinetic energy variations are insignificant when compared to variations in enthalpy or in internal energy. 5. A5.

No moles or energy accumulate in the interface. 6. A6.

The liquid phase, the vapor phase, and the interface operate locally at a state of thermodynamic equilibrium.

It follows from Assumption A1. that molar dynamics for component $j\in\{1,...,c\}$ correspond to

[TABLE]

Under Assumptions A2.-A4., the internal energy can be written using a total energy balance

[TABLE]

and, kinetic energy can be described using a balance on mechanical energy (derived from a momentum balance)

[TABLE]

In equations (4)-(6), $n_{\alpha,j}$, $e_{\alpha}$ and $e_{\alpha,K}$ represent interface sources/sinks for moles (in component $j$) and energy (total and mechanical). Convective inflow rates in (4)-(6), $F_{\alpha,N_{j},\mathrm{in}}$, $F_{\alpha,H,\mathrm{in}}$, $F_{\alpha,V,\mathrm{in}}$, $F_{\alpha,K,\mathrm{in}}$, are assumed as fixed inputs. As there is not a conservation principle for volume (Taylor and Krishna, 1993), we write the liquid volume as a function of the molar holdup (Sandler, 1999)

[TABLE]

where $V_{o}$ is a constant representing the volume of the rigid unit. Volumetric outflow rates are described as a function of the phase flow velocity $v_{\alpha}(\cdot)$, which in turn depends on the kinetic energy of the phase

[TABLE]

where $A_{\alpha,\mathrm{out}}$ stands as a parameter that represents the cross sectional area of the flow line, and $M_{\alpha}=\sum\bar{m}_{j}N_{\alpha,j}$ represents the mass holdup in bulk-phase $\alpha$. Heat flows $Q_{\alpha}$ in the energy balance equation (5) are written as being proportional to differences between the heat exchanger temperature $T_{\alpha,Q}$ and the bulk-phase temperature $T_{\alpha}$

[TABLE]

where $\lambda_{\alpha}$ is a known heat exchange parameter. Temperature in turn is related to the internal energy. The internal energy for phase $\alpha$ can be written as

[TABLE]

where $U_{o,\alpha}(\cdot)=\sum\bar{u}_{\text{o},j}N_{\alpha,j}$ represents the internal energy of the system at a reference temperature $T_{\text{o}}$, and $\mathcal{C}_{\alpha}(\cdot)=\sum N_{\alpha,j}\,\bar{c}_{\alpha,j}$ stands for the total heat capacity for phase $\alpha$. The ideal gas equation is used to write the pressure in the gas phase as a function of the extensive parameters

[TABLE]

Equations (4)-(6) represent conservation principles (Bird et al., 2002). The conservation principles put together with its constitutive equations (7)-(11) stand as a standard representation of an open thermodynamic system. It should be noted here that the description is not yet completed as the interface exchange rates are not included in the constitutive equations (7)-(11). To complete the flash-drum system description, in the next section we describe molar $n_{\alpha,j}$, energetic $e_{\alpha}$ and mechanical $e_{\alpha,K}$ interface exchange rates as function of the intensive variables.

2.2 Interface Exchange Rates

The model presented so far considers the liquid and the gas phases as separated subsystems inside the multiphase unit. As a consequence, the flash-drum is not required to evolve over an equilibrium manifold and inhomogeneities in temperature and chemical potential can arise between phases. The ideas behind the modeling for non-homogeneous systems were first presented by Krishnamurthy and Taylor (1985). Their model relies on describing interfacial temperature, pressure, and compositions solely based on Assumptions A5. and A6. plus a mechanical equilibrium assumption. The non-equilibrium model of Krishnamurthy and Taylor (1985) is unfortunately limited to stationary regimes. Below, we reproduce such description adapted to our open process.

Two boundary layers are assumed to surround the gas-liquid interface as depicted in Figure 3. Liquid and gas coexist in the interface between the layers at local equilibrium. Subindex $i$ is set to describe local interfacial temperature, and molar fractions. In addition, $n$ is defined as the interface molar exchange which represents the total rate at which moles flow from liquid to gas phase. The interface variables are therefore referred as the vector

[TABLE]

Temperature and composition gradients are allowed between the bulk-phases and the interface, see Figure 3. These gradients are known to be the driving forces behind interface flow rates. In the molar balance equation (4), the rate at which chemical component $j$ flows from phase $\alpha$ towards/from the interface is represented by $n_{\alpha,j}$. Neglecting the effects of temperature, pressure gradients, and intercomponent diffusive transport, we can write (Taylor and Krishna, 1993)

[TABLE]

where $y_{j},\;x_{j}$ and $C_{\alpha}$ represent respectively the molar fractions and concentrations in the bulk-phases. The parameter $k_{\alpha,j}$ is assumed to be known and stands for a diffusive transport parameter. The total interface transport rate that appears in (13) corresponds simply to

[TABLE]

The rate at which energy flows from phase $\alpha$ towards/from the interface in Equation (5) is represented by $e_{\alpha}$. The interface energy flow rate can be written as the sum of convective and thermal energy contributions (Taylor and Krishna, 1993; Bird et al., 2002)

[TABLE]

where $\bar{h}_{\alpha,j}$ stands as the partial molar enthalpy of component $j$ in phase $\alpha$, $\lambda_{\alpha,i}$ is a known thermal exchange parameter, $T_{\alpha}$ represents the bulk-phase temperature and $T_{i}$ the interfacial temperature. The rate at which mechanical energy is exchanged between bulk-phase and the interface in the mechanical energy equation (6) is represented as $e_{\alpha,K}$. This term corresponds to (Bird et al., 2002)

[TABLE]

where $\rho_{\,\alpha}$ stands as the density of the bulk-phase $\alpha$. Setting $\bar{m}_{j}$ to be the molar mass of component $j$, we write $m_{\alpha}=\sum\bar{m}_{j}\,n_{\alpha,j}$ for the mass interface transport rate. The term $v_{\alpha,i}$ represents the average velocity of the mass flowing through the boundary layers surrounding the interface. Setting $\text{A}_{i}$ as the interface area we define the interface flow velocity in boundary layer $\alpha$ as

[TABLE]

It should be noted that the interface exchange rates (13)-(15) depend on the $2c+2$ interface variables written in (12). These interface variables are recovered as the solution to an algebraic system of equations built from Assumptions A5. and A6.

2.3 Interface Algebraic System

No accumulation of moles in the interface, Assumption A5., gives $c-1$ equations

[TABLE]

where $n:=n_{l}=n_{g}$ is the total molar interface rate defined in (12). No accumulation of energy in the interface, Assumption A5., leads to one more algebraic restriction

[TABLE]

where $n_{j}:=n_{l,j}=n_{g,j}$ represents the molar interface rate for component $j$ and $\Delta\bar{h}_{\mathrm{vap},j}=\bar{h}_{g,j}-\bar{h}_{l,j}$ represents the partial enthalpy of vaporization for component $j$. To complete the interface description we add $c+2$ interface equilibrium equations

[TABLE]

where $K_{j}(\cdot)$ represents the liquid-vapor composition ratio for component $j$. Even though the liquid-vapor composition ratio is frequently assumed constant, this term is, in general, a nonlinear function of the intensive interface variables, see Equation (57) in B.2.

Balance equations (4)-(6) restricted by the interface equations (16)-(18) sum up to $4c+6$ nonlinear coupled differential algebraic equations that represent a semi-explicit DAE system

[TABLE]

where ${\bf z}$ refers to bulk-phase variables

[TABLE]

The index of a DAE system is the number of times we have to differentiate the constrains $g(\cdot)$ to put aside the algebraic restrictions and write the DAE as an equivalent ordinary differential equations (ODEs) system. A DAE system has index one if and only if the Jacobian $J_{\bf w}(g)$ has full rank (Brenan et al., 1996). It follows from the implicit function theorem that, as $J_{\bf w}(g)$ is non-singular, there exists a unique function $h({\bf z})$ such that substitution of ${\bf w}=h({\bf z})$ satisfies the algebraic restriction in the DAE, i.e., $g({\bf z},h({\bf z}))=0$. Therefore, if a DAE system is of index one we can virtually solve the algebraic part through the mapping ${\bf w}=h({\bf z})$. Substitution of ${\bf w}=h({\bf z})$ into the differential part would transform a DAE model to an equivalent ODE system.

The Jacobian $J_{\bf w}(g)$ for equation (19b) is non singular (see A.1), therefore system (19) is of index one. The irreversible flash-drum DAE model is thus equivalent to an ODE system and it can be characterized using Lyapunov theory. In the next section we put together Gibbs equation (3) with the balance equations in (19a) to compute variations in entropy with respect to time. This forms the basis to calculate the entropy production, a positive definite thermodynamic potential, considered as a Lyapunov function candidate used to characterize the irreversible flash-drum.

3 Entropy Production

In this section, we compute the internal entropy production for the irreversible flash-drum using an entropy balance equation. Being positive definite, the internal entropy production represents a physics-based Lyapunov function candidate which we can use to characterize the dynamics for the nonlinear flash-drum system (19).

Entropy is an extensive non-conserved property (Sandler, 1999). The entropy production rate $\sigma$ for the open flash-drum can be written from an entropy balance as

[TABLE]

where $F_{\alpha,S,\mathrm{in}}$ and $F_{\alpha,S,\mathrm{out}}$ represent convective flow rates of entropy, and $Q_{\alpha}$ stands as a heat source at temperature $T_{\alpha,Q}$. Because entropy is additive over subsystems (Callen, 1985), entropy variations with respect to time can be written using equation (3) as

[TABLE]

where the vector notation

[TABLE]

has been introduced. Note that $\dot{\bf N}_{\alpha}$ and $\dot{U}_{\alpha}$ in (22) represent the molar and energy balances for the flash-drum, equations (4)–(5), and $\dot{V}_{\alpha}$ can be recovered from Equation (7).

Entropy for each phase is a homogeneous of degree one function, Equation (1). Thus, integration of $dS_{\alpha}$, see Equation (2), gives (Callen, 1985, Euler’s theorem)

[TABLE]

Each bulk-phase is perfectly mixed, Assumption A1. Then, multiplying (23) by the inverse of the residence time $1/\tau_{\alpha}=F_{\alpha,V,\mathrm{out}}/V_{\alpha}$ we can write the outflow rate of entropy as

[TABLE]

Substitution of (22) and (24) into (21) gives the entropy production for the irreversible flash-drum as

[TABLE]

Each term in Equation (25) can be identified as a source of entropy on a phenomenological basis:

• •

The first terms,

[TABLE]

represent the entropy produced as heat $Q_{\alpha}$ is exchanged between phase $\alpha$ and the external heat sources.

• •

The terms related with convective flow rates,

[TABLE]

represent the entropy produced as the inflows are mixed with the bulk-phases. Note that if inflows are at the same temperature and chemical potential (composition) as the bulk-phases, these terms are equal to zero.

• •

The last two terms,

[TABLE]

represent the entropy produced as energy and moles are exchanged through the interface when phases are not homogeneous, i.e., far from thermodynamic equilibrium.

Note that when the flash-drum is completely isolated from the environment we have $F_{\alpha,H,\mathrm{in}}=0$, ${\bf F}_{\alpha,N,\mathrm{in}}=0$ and $Q_{\alpha}=0$. Then, the entropy production is reduced to the internal entropy production

[TABLE]

Equation (26) can be written as the sum of products between generalized flows and driving forces

[TABLE]

where the flows, see equations (13), (16) and (14),(17), correspond to

[TABLE]

and the driving forces

[TABLE]

come as a consequence of inhomogeneities between liquid and gas phases. Even though the flows (27) and the driving forces (28) are not written using Onsager’s relations (Prigogine, 1968), it can be easily verified that both (27) and (28) vanish at thermodynamic equilibrium. Moreover, numerical evidence shows that $\sigma_{i}$ is positive definite and decreases with respect to time as the flash-drum reaches an equilibrium state in a system with fixed inflows, as shown in the next section.

4 Numerical Simulations

In this section, a non-ideal methanol-water mixture is simulated to illustrate the properties of the model presented in Section 2 and Section 3. First, a change of coordinates is used to transform the extensive-intensive description (19) to an equivalent DAE system improved for numerical integration. Then, a stability analysis based on Lyapunov first method is briefly discussed. Numerical trajectories demonstrate that the internal entropy production can be considered as a Lyapunov function candidate to assess the stability for the nonlinear irreversible flash-drum model.

4.1 Change of Coordinates

Extensive variables in (19a) are considerably larger in magnitude than the intensive variables in the interface description (19b). Such differences are known to cause the Jacobian matrix of the system to be ill-conditioned (Ritschel et al., 2018). This leads to precision problems during the numerical integration of a DAE or an ODE system. To avoid scale differences between the bulk-phase and the interface models, we use standard definitions for molar fractions and molar concentration together with constitutive equations (7)-(10) to write a change of coordinates

[TABLE]

given by

[TABLE]

Here $N_{\alpha}$ and $M_{\alpha}$ represent the total molar and mass holdups for phase $\alpha$, $U_{\alpha,o}$ is a reference state at temperature $T_{o}$, $\mathcal{C}_{\alpha}$ is the heat capacity for phase $\alpha$, $V_{o}$ is the volume of the flash-drum, and $\bar{v}=\sum\bar{\text{v}}_{j}x_{j}$ represents the liquid molar volume.

As the change of coordinates is bijective (A.2), time differentiation of (30) leads to a description of the flash-drum equivalent to (19) with the advantage that the bulk-phase and the interface variables have the same order of magnitude. Taking a time derivative on both sides of (30) allows us to rewrite the DAE (19) as

[TABLE]

where

[TABLE]

represents the new bulk-phase variables, and

[TABLE]

denotes again the interface variables. The matrix M in (31) is diagonal of dimension $2c+4$

[TABLE]

where $1_{c-1}$ stands for a row vector of dimension $c-1$ with ones as elements. The entries in the vector function $h$ in (31) are given by

[TABLE]

where sub-index “$\mathrm{in}$” refers to inflows, $\bar{\mathcal{C}}_{\alpha}$ represents molar heat capacity for the flow $\alpha$, $m=\sum\bar{m}_{j}n_{j}$ holds for the interface mass exchange rate, $n_{l,V}=\sum\bar{\text{v}}_{j}n_{j}$ represents the volumetric change rate at which the liquid phase evaporates, and $P_{g}n_{g,V}=n\,R\,T_{g}$. The algebraic restriction $g$ in (31), defined before in (19), is rewritten here for the sake of clarity

[TABLE]

Equation (31) represents again an index one DAE system describing the irreversible flash-drum. The system has $2c+4$ bulk-phase variables and $2c+2$ interface variables. Conservation principles applied to each bulk-phase give $2c+4$ differential equations, while the interface is described by $2c+2$ algebraic equations. In contrast with (19), all the variables in (31) have the same order of magnitude, reducing the possibility of numerical integration problems.

4.2 Local Stability

In order to study the dynamic properties of the proposed model, we briefly review local stability results for linear index one DAE systems. A complete view on stability for DAE systems can be found in the work by Yang et al. (2013).

Let $({\boldsymbol{z}}^{\star},{\bf w}^{\star})$ be an equilibrium point of (19), i.e.,

[TABLE]

As the Jacobian $J_{\bf w}(g)$ is full ranked (see A.1), we can write a linearized version of (31) around as $({\boldsymbol{z}}^{\star},{\bf w}^{\star})$ as

[TABLE]

where the matrix $\Theta$ corresponds to the linearization of the right hand side terms in (31a) at the equilibrium point

[TABLE]

Lyapunov’s first method states that if the spectrum of the matrix (34) is contained in the left half-complex plane, then (31) is locally asymptotically stable in a neighborhood of the stationary state $({\boldsymbol{z}}^{\star},{\bf w}^{\star})$.

4.3 Case Study

A non-ideal methanol-water mixture is considered to illustrate the proposed model and its analysis. With two components, the DAE system (31) is given by 14 equations and 14 variables. Fixing the inflows to be at thermodynamic equilibrium $(T^{\star},\,P^{\star},\,y^{\star}_{1},\,y^{\star}_{2},\,x^{\star}_{1},\,x^{\star}_{2})$, the stationary state for (31) corresponds to

[TABLE]

where $C^{\star}_{g}$ is determined from Equation (11a), $V^{\star}_{l}$ is the volume of the liquid phase at equilibrium, and $F_{\alpha,V}^{\star}$ represents the stationary state volumetric inflow rates. To determine numerical values for the steady state, inflow properties are assumed to be at thermodynamic equilibrium at $T^{\star}=351.24\mathrm{K}$ (78.09${}^{o}\mathrm{C}$) and $P^{\star}=101.3\mathrm{kPa}$ (0.9998 $\mathrm{atm}$). Liquid-vapor equilibrium is calculated using an Antoine–Margules thermodynamic model (see B.2). In addition, inflows are fixed at $1\;\mathrm{m}^{3}/s$, and the liquid phase is set to occupy 10% of the total volume $V_{T}=1\;\mathrm{m}^{3}$. At these conditions, the stationary state $({\boldsymbol{z}}^{\star},{\bf w}^{\star})$ takes the values given in Table 1.

4.3.1 Stability Analysis

Below we assess local stability when the system (31) is linearized through (33) at the stationary state given in Table 1.

1. S1

Scenario 1 (non-isobaric operation regime). At equilibrium, the linearized system (33) has rank 8. This system is unstable as the spectrum of the system has two positive eigenvalues $\lambda_{1}\approx 7.7691\times 10^{-3}$, and $\lambda_{2}\approx 1\times 10^{-12}$. The remaining eigenvalues are contained in the left half plane, between $\lambda_{3}\approx-1$ and $\lambda_{8}\approx-1.2\times 10^{5}$. 2. S2

Scenario 2 (isobaric operation regime). We calculate the stability properties of the system after a perfect pressure controller is included in (31) making pressure constant at $P_{l}=P_{g}=P^{\star}$. It must be noted here that the introduction of the perfect pressure controller changes the dynamic description. Gas concentration trajectories are described under an isobaric regime by the algebraic restriction

[TABLE]

where functions $h_{2c-1}$ and $h_{2c+3}$ represent the right hand terms for the temperature and concentration derivatives in (31a):

[TABLE]

Equation (35) can be solved to write $C_{g}$ as a state function. Then, the dimension of the linearized flash-drum is reduced as $C_{g}$ is an equation of state for the isobaric system and not a state variable as in the non-isobaric model (31). Under the isobaric restriction (35) the linearized system (33) has rank 7.

One eigenvalue for (33) restricted by the isobaric equation (35) seems to be positive, $\lambda_{1}\approx 1\times 10^{-13}$. The rest of the spectrum is contained in the negative half-line between $\lambda_{2}\approx-1$ and $\lambda_{7}\approx-1.2\times 10^{5}$. The positive eigenvalue is so close to zero that we cannot draw a sound conclusion regarding stability as numerical error can be the underlying cause for the positivity of $\lambda_{1}$.

4.3.2 Numerical Simulations for Isobaric Operation

To further investigate the stability of the isobaric regime, three dynamic simulations for Scenario 2 are now discussed. For the first scenario, the system is initialized at stationary state and the inflow liquid temperature is reduced, pushing the system far from the thermodynamic equilibrium state. In the second simulation, the system is disturbed from the stationary state. The disturbance is then removed at time $t=2\mathrm{s}$ and the system goes back to the equilibrium state. For the last simulation we show how the internal entropy decreases as the system goes back to thermodynamic equilibrium starting from non-equilibrium initial states.

Numerical simulation 1 (Figure 4). The inflow liquid temperature is disturbed through a ramp disturbance for $1\leq t<2$. The liquid inflow reaches a temperature $T_{l,\mathrm{in}}=0.95\times T^{\star}$ for $t\geq 2$ and the liquid-vapor system is forced to operate far from thermodynamic equilibrium. As the system remains far from thermodynamic equilibrium, mass and energy flow from the gas to the liquid phase. The gas phase condensates completely around $t\approx 460\mathrm{s}$. For better appreciation of the dynamic behavior Figure 4 is presented using a log scale for time.

Numerical simulation 2 (Figure 5 and Figure 6). A ramp disturbance is introduced in the liquid inflow temperature for $1\leq t<2$. In contrast with the numerical simulation 1, the disturbance is removed at $t=2\mathrm{s}$ and the inflow temperature is set back its nominal value $T_{l,\mathrm{in}}=T^{\star}$, see Figure 5. This scenario is studied through two different disturbances. It can be seen in Figure 6 that inhomogeneities in temperature and composition appear between phases as a consequence of the disturbance. Then, transfer processes redistribute the mass and the energy in the system as the flash-drum goes back to the equilibrium state.

Numerical simulation 3 (Figure 7). Results from the Numerical simulation 2 point towards the equilibrium state being a stable steady state (Figure 6), despite the presence of a zero eigenvalue in the linearized system (see Section 4.3.1). To extend the analysis, the last simulation here presented tests the internal entropy production (26) as a Lyapunov function candidate for the irreversible flash-drum111Dynamic trajectories on the internal entropy production (26) are computed using the ideal gas model and the fugacity to calculate the chemical potentials (B.1).. In this scenario, the liquid-vapor system starts from an initial condition far from thermodynamic equilibrium while the inflows remain constant at the conditions in Table 1. It can be seen that the internal entropy production, Equation (26), behaves as a Lyapunov function as the irreversible flash-drum goes back to the stationary state, see Figure 7. This suggests again that the thermodynamic equilibrium state is a stable stationary state for the irreversible flash-drum.

5 Conclusions and Future Work

In this article, modeling aspects of dynamic flash-drum systems are explored using a non-equilibrium physics-based model. The description here presented considers transport phenomena as well as conservation principles to write the dynamics of a multiphase system as a nonlinear DAE system of index one. The proposed dynamic model describes the evolution of liquid and vapor phases as separated sub-systems interconnected through an interface. The introduction of the interface exchange rates in the model can even predict the collapse of one phase for systems that operate consistently far from equilibrium. Moreover, the model presented here describes how entropy is produced as a consequence of external (mass/energy exchanges between the system and the environment) and internal phenomena (mass/energy exchanges between phases).

Numerical evidence shows that the linearized irreversible flash-drum system has positive eigenvalues for a non-isobaric operation regime. For the isobaric case, trajectories appear stable in simulations as the internal entropy production appears to be a Lyapunov function candidate for the system. The numerical results pinpoint the need to perform a deeper analysis regarding the stability for the nonlinear liquid-vapor irreversible system. In future research, a nonlinear passivity-based stability analysis approach for the DAE model along the lines of García-Sandoval et al. (2015) will be considered to get an input-output perspective on the analysis and control problem for multiphase chemical systems.

Acknowledgments

The research presented in this paper is supported by the Mexican Council for Science and Technology (Grant 410828), and the Institute of Information and Communication Technologies, Electronics and Applied Mathematics (ICTEAM) at Université catholique de Louvain, Belgium.

Appendix A Jacobians

A.1 Jacobian for the algebraic system g

The Jacobian of $g$, see Equation (19), with respect to

[TABLE]

corresponds to the full rank222The Jacobian matrices and the respective ranks in A are calculated using the open source computer algebra system Maxima version:5.32.1. $2c+2$ square matrix

[TABLE]

The term $\text{k}_{\alpha}$ in (45) represents a $(c-1)\times c$ diagonal sub-matrix with an additional column of zeros

[TABLE]

The terms ${\text{y}}_{c}$ and ${\text{x}}_{c}$ stand for a column vectors

[TABLE]

When subindex $c-1$ is used to write ${\text{y}}_{c-1}$, and ${\text{x}}_{c-1}$, the composition vector contains only the first $c-1$ molar fractions. The term $\Delta\text{h}_{\alpha}$ represents a row vector

[TABLE]

and the vector $\text{h}_{\alpha}$ contains the enthalpies for the components in the mixture

[TABLE]

The symbol $\text{x}_{i}$ stands for a diagonal matrix with interface molar compositions

[TABLE]

and $J_{\beta}(K)$ holds for the Jacobians of the equilibrium ratio $K$ with respect to $\beta\in\{y_{1,i}\ldots y_{c,i}\;\;x_{1,i}\ldots x_{c,i}\;\;T_{i}\;\;n\}$. Zeros and ones in (45) represent zero vectors, one vectors and scalars when dimension is not specified.

A.2 Jacobian for the change of coordinates

To show that that Equation (30) is bijective, it is enough to demonstrate that the Jacobian matrix of the mapping is non-singular. To do so, we rewrite the change of coordinates as

[TABLE]

where

[TABLE]

Then, the Jacobian of $F$ with respect to

[TABLE]

corresponds to the full rank sparse matrix

[TABLE]

where, $\Theta_{g}(\cdot)$ stands for a $(c-1)\times c$ matrix

[TABLE]

and $\Theta_{l}(\cdot)$ holds for an equivalent matrix in terms of liquid molar fractions. The molar energy, the liquid molar volume, and the molar mass vectors in (51) correspond to

[TABLE]

Finally, the zeros and ones in (51) represent zero matrices/vectors, one vectors and scalars when dimension is not specified.

Appendix B Chemical potential and thermodynamic equilibrium

B.1 Chemical potential for gas and liquid mixtures

Consider a thermodynamic system formed by an ideal gas with internal energy $U$, and $N$ moles occupying a volume $V$. The entropy for that system satisfies Callen (1985)

[TABLE]

where the reference state $\bar{\text{s}}_{\text{o}}$ corresponds to the molar entropy for a system at temperature $T_{\text{o}}$ with molar volume $\bar{\text{v}}_{\text{o}}$, and the dimensionless heat capacity $\mathtt{c}_{\text{v}}=c_{v}/R$ is considered a constant parameter. Values for $\mathtt{c}_{\text{v}}$ are reported in Table 2.

For a mixture of ideal gases occupying a volume V at temperature $T$, entropy is the sum of entropies that each component would have if it alone were to occupy the volume V at temperature $T$ (Callen, 1985, §3.4 – Gibbs theorem). Then entropy for an ideal mixture of gases can be written as

[TABLE]

where the sum is taken over all the components in the mixture, $N=\sum N_{j}$ stands for the total molar holdup of the system, and $\mathcal{C}=\sum\mathtt{c}_{\text{v},j}N_{j}$ represents the total heat capacity of the system. Note that Equation (52) is the entropy fundamental equation defined by the thermodynamics formal structure (1) for an ideal gas mixture

[TABLE]

It follows that chemical potential for component $j$ inside an ideal gas mixture can be calculated as the derivative of (52) with respect to mole numbers

[TABLE]

Despite Equation (53a) being only valid for ideal gas mixtures, it is a common practice to try to preserve this form as far as possible when describing non-ideal systems Prigogine (1968).

The chemical potential for a component $j$ inside a non-ideal mixture can be written as Prigogine (1968)

[TABLE]

where $\mu^{\mathrm{im}}$ represents the chemical potential of component $j$ inside an ideal mixture. To measure deviations from ideal behavior, excess in chemical potential (referred by some authors as excess on partial molar Gibbs potential $\bar{g}_{i}^{\mathrm{ex}}$) is defined as

[TABLE]

and thus the activity coefficient $\gamma_{j}$ is satisfies

[TABLE]

Setting the chemical potential for $j$ inside the ideal system $\mu_{j}^{\mathrm{im}}$ to be represented by (53) we can write the chemical potential for component $j$ inside a liquid mixture as

[TABLE]

This description for the chemical potential has a clear physical interpretation. The chemical potential $\mu_{j}$ is the chemical potential of $j$ as an ideal system (an ideal gas in this case), plus mixing effects (second term), plus a correction term that represents deviations from ideal behavior $\mu^{\mathrm{ex}}$. In the following section, we write the excess in chemical potential $\mu_{j}^{\mathrm{ex}}$ using a liquid-vapor equilibrium model.

B.2 Margules-Antoine equilibrium model

A liquid-vapor system with two components at temperature $T$ and pressure $P$ is said to be at thermodynamic equilibrium when the chemical potentials in liquid an gas phases are equal

[TABLE]

As chemical potential is not measurable, it is common practice to rewrite (55) as an equivalent equality between fugacities

[TABLE]

where $\bar{f}_{j}^{\,g}=y_{j}P$, provided that $P$ is close to atmospheric pressure, stands as the gas fugacity for component $j$ in the gas mixture, and

[TABLE]

Parameters $A_{j},\,B_{j},$ and $C_{j}$, in the previous equation stand as constants for Antoine’s Equation. The excess in chemical potential can be written as a function of a polynomial in the liquid composition $Q(x)$ (Taylor and Krishna, 1993)

[TABLE]

Given thermodynamic parameters $A_{12}$ and $A_{21}$, the liquid composition polynomial $Q$ can be written using Margules thermodynamic model as

[TABLE]

Algebraic rearrangement of equation (56) leads to the nonlinear liquid-vapor equilibrium equation

[TABLE]

where

[TABLE]