The influence of interfacial transfer and film coupling in the modeling of distillation columns to separate nitrogen and oxygen mixtures

https://doi.org/10.1016/j.cesx.2020.100076 2590-1400/ 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). ⇑ Corresponding author at: Universidad de Buenos Aires, Facultad de Ingeniería, Departamento de Química, Paseo Colón 850, C1063ACV Buenos Aires, Argentina. E-mail address: dkingston@fi.uba.ar (D. Kingston). Diego Kingston a,d,⇑, Øivind Wilhelmsen , Signe Kjelstrup d


Introduction
Evaporation and condensation are examples of first order phase transitions that occur in process equipment such as distillation columns (Krishnamurthy and Taylor, 1985), chemical reactors (Wilhelmsen et al., 2018) and multiphase heat exchangers (Wilhelmsen et al., 2014), as well as in natural processes such as cloud formation (Kulmala et al., 2004) and evaporation from oceans (Neale et al., 2010). Accurate prediction of heat and mass transfer rates is a prerequisite for reliable modeling of these processes. A simple description of mass transfer is provided by Fick's law, which dictates that the diffusion flux goes in the opposite direction of the concentration gradient (Taylor and Krishna, 1992). Fourier's law is an analogous description of heat transfer that prescribes that the heat flux goes in the opposite direction of the temperature gradient . Supplementary to these laws, Ludwig (and later Soret) showed in 1856 that a temperature gradient can induce a mass flux (Köhler and Morozov, 2016;Ludwig, 1856). This is called thermodiffusion or the Ludwig-Soret effect. Conversely, a concentration gradient can induce a heat flux (the Dufour effect) . These phenomena are linked and referred to as coupled transport processes (Kjelstrup and Bedeaux, 2008).
Non-equilibrium thermodynamics (NET) provides a consistent and general framework to deal with coupling : The fluxes are functions of all thermodynamic driving forces, J k ¼ P L ki X i , with J k being the fluxes, X i , the forces and L ki , the phenomenological coefficients. Here, L kk are the main coefficients (those linking a flux with its main driving force), while the L ki are coupling coefficients that obey Onsager's reciprocal relations, L ki ¼ L ik . The latter are related to the Soret and Dufour coefficients in the case of coupling between heat and mass transfer. Other examples where coupled transport processes are important are thermoelectric generators Savani et al., 2017) and modules to extract power from mixing pure water with sea water by use of reverse electrodialysis (Zlotorowicz et al., 2017).
In bulk-systems, the coupling coefficients are usually less important than the main coefficients. Across interfaces, however (Kjelstrup and Bedeaux, 2008), coupling is of higher relevance due to the latent heat associated with the phase transition (Wilhelmsen et al., 2014). It has for long been known that the interface poses a resistance to heat and mass transport, as noticed by Kapitza in the observation of a temperature "jump" across a https://doi.org/10.1016/j.cesx.2020.100076 2590-1400/Ó 2020 The Author(s). Published by Elsevier Ltd. This is an open access article under the CC BY license (http://creativecommons.org/licenses/by/4.0/). liquid-solid interface (Wilhelmsen et al., 2014). The resistances have been quantified by kinetic gas theory, statistical rate theory, non-equilibrium molecular dynamics simulations, density functional theory and experiments (Glavatskiy and Bedeaux, 2010;Klink et al., 2015;Lamanna et al., 2018;Wilhelmsen et al., 2014). Kinetic gas theory describes the behavior of the interface in terms of velocity distribution functions, and has resulted in the wellknown Hertz-Knudsen equation and its extensions which predict the net mass flux during evaporation and condensation (Bond and Struchtrup, 2004). These formulas include the evaporation (or condensation) coefficient as an empirical parameter. When calculated from experiments with nominally the same experimental conditions, the evaporation coefficients differ widely from one another (Tsuruta and Nagayama, 2002). Statistical rate theory was developed to predict the net mass flux during evaporation (Ward and Fang, 1999). The evaporation rate is here expressed in terms of an evaporation probability, defined through the entropy change of a molecule going from the liquid phase to the vapor phase. Recent discussions have pointed out that statistical rate theory does not provide an expression for the interfacial energy flux, which plays a leading role in determining the interface temperature-jump as discussed in the literature (Badam et al., 2007). For most systems apart from the truncated and shifted Lennard-Jones fluid (Wilhelmsen et al., 2015) and water (Wilhelmsen et al., 2016), the interface transfer coefficients required in NET are unknown parameters.
The coupling of transport processes, as well as the resistance posed by the interface, are usually neglected in the modeling of both natural processes and industrial processes (Wilhelmsen et al., 2016). The validity and consequences of these assumptions are mostly unknown. There is a need for further insight on when these phenomena should be accounted for. Studies by van der Ham et al. (2010) and Mendoza and Kjelstrup (2011) have shown that the interfacial transport has an important effect in the prediction of the mass and measurable heat fluxes at the level of the interfacial region. However, the influence of these phenomena on the overall performance of a process unit is still unresolved. In this work, we present a rigorous assessment of the relevance of these phenomena in the modeling of a distillation column for air separation, which we consider to be a binary mixture of nitrogen and

Latin symbols a
Interfacial area between liquid and vapor per unit volume of packing m 2 Á m À3 À Á  oxygen. In more realistic modeling of air separation, argon should also be included in the mixture. This example has been chosen because it is an important process in its own right (van der Ham, 2011) and because its overall performance is highly influenced by the interchange of heat and mass between the liquid and vapor phases. Distillation is one of the most widely used and energy consuming separation operations (Kingston et al., 2020). This has motivated numerous studies to enhance its performance, both with respect to the energy and mass transfer efficiency, and to develop models for simulation, both at steady state and in dynamic operation. For example, Meng et al. (2019) studied SiC ceramic foam as packing in distillation columns. They performed pilot plant experiments for two binary mixtures (nitrogen-oxygen and oxygen-argon) and found that SiC foam provides higher operation flexibility and significantly enhances mass transfer efficiency. Rizk et al. (2012) performed an exergy analysis of a packed heat-integrated cryogenic air separation unit using a rate-based model and showed that its energy efficiency is 23% higher than that of a conventional double column. This was achieved by operating the high-pressure column at lower pressures, which resulted in lower duty in the condenser and a lower compression power demand. Their model had some limitations, as they considered only diffusive transport. Chang and Liu (2014) used a more rigorous rate-based model to simulate a heat-integrated distillation column (HIDiC) to separate nitrogen, oxygen and argon, estimating the total capital and operating costs for this column and showing that the vapor Murphree efficiencies were below 0.7. Wang et al. (2019) evaluated different mass transfer correlations for cryogenic air separation by performing simulations with a rate-based model and comparing the results with experimental data from a structured packing column with a capacity of 17,000 Nm 3 Á h À1 . They found that, for 3 different groups to predict the mass transfer coefficients, the differences in nitrogen mole fraction were below 2.7%; the authors concluded that the use of rate-based models with these correlations is appropriate to simulate the steady state characteristics of cryogenic air separation. Kingston et al. (2020) performed an entropy production minimization in diabatic distillation columns to separate a binary mixture of nitrogen and oxygen, using rate-based and equilibrium-stage models. They searched for the optimal heating/ cooling strategy and found that the entropy production can be reduced by nearly 50%, while decreasing the total heating and cooling duties by 30 and 50%, respectively. They also reported that, for the equilibrium-stage model, there is a trade-off between lower entropy production and higher heating duties. Zhou et al. (2019) developed a rate-based equation-oriented parallel column model to simulate divided wall columns, which is compliant with any CAPE-OPEN process flowsheet simulation package. They validated their model with experimental data from pilot plant measurements for alcohol and hydrocarbon systems, finding very good agreement in the temperature profiles and in the prediction of the key component's mole fraction. Lashkarboluki and Mallah (2019) used a rate-based model to study the dynamic behavior of a distillation column for 18 O and compared their results with experimental data from a semi-industrial column, with an estimated error of 5%. State-of-the-art, mathematical models for distillation columns rely on the assumption that the interface is in thermodynamic equilibrium (Taylor and Krishna, 1992;Taylor et al., 2003). Lamanna et al. (2018) found, by use of Raman measurements, that models assuming equilibrium at the interface failed to provide a good estimate of the temperature of n-hexane droplets. This indicates that the commonly used models have potential for improve-ment by considering also the interfacial resistances. Furthermore, the coupling between heat and mass transfer in the films adjacent to the interface is also usually neglected, e.g. in the wellestablished model for distillation columns by Taylor and Krishna (1992) and Taylor et al. (1994). They state that for classical unit operations (e.g. distillation), the temperature gradients are rarely high or sustained long enough to make thermal contributions to mass fluxes important. In this work, we shall quantify these effects in further detail for the separation of nitrogen and oxygen mixtures. Kjelstrup and de Koeijer (2003) and van der Ham et al. (2010) extended the film model to describe transport between vapor and liquid with NET (Taylor and Krishna, 1992). They modeled the interfacial region as a series of connected control volumes and obtained expressions for the overall resistivities. van der Ham et al. (2010) studied the effect of these phenomena for nitrogen and oxygen binary mixtures, Mendoza and Kjelstrup (2011), for ethanol and water and Keulen et al. (2017) modeled membrane distillation for the production of drinking water. A detailed discussion of the relevance of these transport mechanisms for the overall performance of the process unit is still missing. This could be important in the field of process design. Using a distillation column for nitrogen and oxygen separation as example, we shall give insight into when and how interfacial resistance and coupled transport in the adjacent films become important.
The article will be organized as follows. In Section 2, we present the governing equations and assumptions for the distillation column model, the detailed model for transport across the vapor-liquid interface and the adjacent vapor and liquid films. Next, in Section 3 we elaborate on the numerical methods used in the solution of the interface conditions and the distillation column. In Section 4 we define cases to evaluate the relative importance of the different transport mechanisms. Results are presented in Section 5 and concluding remarks are given in Section 6.

Description of the system and model
In this section, we present a mathematical model of a distillation column which separates a binary mixture of nitrogen and oxygen (Section 2.1 and a detailed model for transport across the vapor-liquid interface and the adjacent gas and liquid film phases (Section 2.2). The purpose of these models is to enable an assessment of the relevance of transport across the vapor-liquid interface and coupled transport in the films on the overall performance of the process unit. Boundary conditions for the mathematical models are presented in Section 2.3 and a consistency check of the model based on the entropy production is elaborated in Section (2.5). The Peng-Robinson equation of state was employed to calculate the thermodynamic properties of the mixture. The performance was evaluated by comparing to experimental data for the oxygen-nitrogen mixture (van der Ham, 2011), which showed that there was maximum 5% error in the calculation of molar enthalpy and around 1% for bubble pressure and vapor composition in the temperature, pressure and composition ranges considered. The components' critical properties, acentric factors and ideal gas heat capacities were taken from Perry and Green (2008).

The distillation column model
In the following, we present a mathematical model of a continuous distillation column for separation of a binary mixture of nitrogen-(1) and oxygen-(2), following Johannessen (2004) and Taylor and Krishna (1992). The column includes a total condenser and a partial reboiler (Seader et al., 2011). A sketch of the process unit is shown in Fig. 1, while a more detailed picture of the mass transfer process is given in Fig. 2. The model relies on the following assumptions: 1. The column operates at steady state. 2. The components' mole fraction, velocity and temperature of each phase are uniform in a cross-section. 3. There is no entrainment of liquid in the vapor phase, nor is there vapor trapped in the liquid. 4. There is no back-mixing in the column. 5. There are no chemical reactions. 6. The pressure drop across the section and auxiliary equipment (reboiler and condenser) is negligible.
7. The reboiler behaves as an ideal equilibrium-stage. This means that the outgoing vapor and liquid streams are in thermodynamic equilibrium with respect to each other.
The governing equations for the rectification and stripping sections were given by Johannessen (2004), The conservation of mass for each component in the vapor and liquid phases is given by Eqs. (1) and (2), and the energy balances for the vapor and liquid are given in Eqs. (3) and (4) To solve this system of differential equations, it is necessary to provide an additional model to compute the fluxes from the vapor to the liquid phase.

Model for transport of heat and mass across the vapor-liquid interface and the adjacent film phases
Since the fluids close to the interface are nearly stagnant, part of the resistance to heat and mass transfer originates in thin boundary-layers/films located on each side of the interface. To obtain heat and mass fluxes between the liquid and vapor phases,  we will employ a commonly used representation that considers the overall resistance to come from three resistances in series; a liquid film, the interface and a vapor film (van der Ham et al., 2010). An illustration of this representation is shown in Fig. 1. Furthermore, the following assumptions are used in developing the local transport model, 1) The system is at steady state, 2) The pressure is constant across the film-interface-film system. With these assumptions, the molar and energy fluxes are constant (Taylor and Krishna, 1992). The energy conservation across the system relates the measurable heat fluxes at points u 0 and u 1 by means of At a given location in the distillation column, this equation enables us to calculate the measurable heat flux at any position across the film-interface-film system, u, if the value of the other measurable heat flux is known, or alternatively, if the energy flux, pressure and local temperature are known. We shall next present a model based on NET for describing transport of heat and mass across the system. If the coupling coefficients for heat and mass transfer in the bulk phases are zero and the resistivities at the interface are neglected in this model, we recover Taylor and Krishna's model (Taylor and Krishna, 1992). Therefore, the model presented in this work generalizes previous treatments and enables an assessment of effects that were not previously taken into account.

Transport across the vapor and liquid film phases
To describe transport across the films, we use the force-flux relations from NET van der Ham et al., 2010), where r a ij is the resistivity linking force i with flux j in phase a; l a k is the chemical potential of component k in phase a and the subscript T indicates that the derivative is taken at constant temperature. The resistivities can be cast into a symmetric matrix due to Onsager's reciprocal relations (Kjelstrup and Bedeaux, 2008). This matrix is singular, because the thermodynamic forces are not independent due to the Gibbs-Duhem equation (Kjelstrup and Bedeaux, 2008). Since we work with two components, nitrogen and oxygen, the resulting force-flux relations after imposing the Gibbs-Duhem relation become: where R is the universal gas constant, z a 1 , the mole fraction of nitrogen in phase a (subscript 2 refers to oxygen), v is a dimensionless position, d a is the film thickness in phase a and C a 1 is the thermodynamic factor given by, where/ a 1 is the fugacity coefficient of nitrogen in phase a. These are the governing equations for the film phases. The resistivities in the bulk phases can be calculated from typical transport coefficients, such as the thermal conductivity, k a , the Maxwell-Stefan diffusivity, D a , (which is convenient to use when both components are taken into the description as here) and the heat of transfer, q I;a i . The expressions relating the resistivities and these coefficients were presented in the work by van der Ham et al. (2010): with c a being the molar density of phase a. We used the same constant values for the Maxwell-Stefan diffusivities, heat of transfers and thermal conductivities as van der Ham et al. (2010).

Transport across the vapor-liquid interface
The vapor-liquid interface is also a resistance to heat and mass transfer (Kjelstrup and Bedeaux, 2008). The force-flux relations of the interface are (Kjelstrup and Bedeaux, 2008): 1 where I indicates that the properties are those of the interface and the chemical potentials are evaluated at the temperature of the liquid side of the interface. Here, the R ij are the interface transfer coefficients, which obey the Onsager reciprocal relations (Kjelstrup and Bedeaux, 2008). Since no experimental values or simulation results for the interface transfer coefficients are available for the oxygennitrogen mixture, we used the values given by van der Ham et al. (2010), which are based on kinetic theory. The uncertainty of the predictions from kinetic gas theory is unknown; this will be further discussed in Section 5.

Boundary conditions
2.3.1. Interface problem At each location in the distillation column, the bulk temperature, T a z ð Þ, and the composition, z a 1 z ð Þ, of each phase result from solving the differential equations presented in Section 2.1. This provides the necessary boundary conditions at Locations 1 and 4 in Fig. 1. In particular, we have that The

Distillation column
The feed introduces a discontinuity in the temperature and molar flow profiles. Therefore, the column model must be divided into two different regions (rectification and stripping sections). The appropriate boundary conditions must be provided to connect the regions where the feed is introduced. The feed location has been fixed to z ¼ 0, with the rectification section spanning from z ¼ 0 þ to z ¼ ' R , and the stripping section spanning from z ¼ À' S to z ¼ 0 À . The solutions at z ¼ 0 þ and z ¼ 0 À must satisfy the following mass and energy balances, In the above equations, F k is the component molar flow in the feed, / is the vapor fraction of the feed stream and h a is the molar enthalpy of phase a, evaluated at the corresponding temperature and composition. At the top of the rectification section, there is a condenser. The temperature of the reflux is assumed to correspond to the bubble temperature of the mixture. If the reflux ratio, r, is specified, we have that where T B is the bubble temperature at the condenser pressure and composition of the reflux. A bubble point calculation must be performed when solving these boundary conditions. We consider the reboiler to be an ideal equilibrium-stage. Therefore, the vapor's temperature corresponds to the dew point value, T D . This determines the composition of the liquid, x I k , in equilibrium with the vapor. If the bottoms flow rate, B, is given, the following boundary conditions must be satisfied:

Application of this model to tray columns
The model we have presented applies to columns with random or structured packing. Taylor and Krishna (1992) indicate that packed columns can be simulated as tray columns by taking a height of packing as one stage/tray. Therefore, this model may be used, with some modifications, to tray columns. The changes apply to the calculation of heat and mass transfer rates and the equilibrium conditions at the interface, which are replaced by the differential equations for the vapor and liquid films, Eqs. (9) and (10) (Taylor and Krishna, 1992) for details. This differential-algebraic system of equations can be converted into a nonlinear algebraic system by discretization of the differential equations using a coarse grid and solved by usual techniques for tray columns, for example, those presented in Taylor and Krishna (1992).

The entropy production
The entropy balance facilitates a consistency check of the model. The entropy balances in the films are, where J S is the entropy flux and r is the local entropy production (per unit volume). After integrating the above expression across each film phase and including the contribution from the interface, we obtain with r 0 I being the local entropy production at the interface (per unit area). If we replace the entropy flux by its contributions, we obtain where s a k is the partial molar entropy of component k in phase a.
This equation indicates that the difference in the entropy fluxes into and out of the system is equal to the entropy production per unit area, r 0 , The local entropy production can also be calculated from the forces and fluxes, An analogous equation can be formulated for the entropy production at the interface (Kjelstrup and Bedeaux, 2008). If the model is consistent and well implemented, the two different ways of computing the entropy production give the same result. The local entropy production (per unit area) obtained at each point can be integrated over the interfacial area in the column to give the total entropy production in the section, where the integral is performed over the total interfacial area in the section, A Sec . It should give the same value as the macroscopic entropy balance in the section (Kingston et al., 2020;Kjelstrup et al., 2017).

The numerical procedures
3.1. Solution method for the film-interface-film system The boundary value problem given by Eqs. (7) and (8) was solved using MATLAB's solver bvp4c with three unknown parameters. The relative tolerance was set to 10 À6 . The molar fluxes, J 1 and J 2 , and the measurable heat flux in the vapor, J 0;V q were defined as unknown parameters. The initial solution was constant in each film, using the values corresponding to the bulk conditions. This was the procedure used to solve the film-interface-film system in Section 5.2.

Solution strategy for the complete distillation column
The boundary value problem (BVP) given by Eqs. (1)-(4) was solved using MATLAB's solver bvp4c, which exploits a collocation method, with 10 À3 as relative tolerance. At any location in the distillation column, the molar and heat fluxes must be computed. As solving a boundary value problem at each location in the distillation column is computationally very expensive, we solved the coupled problem involving Eqs. (7), (8), (18)-(20) in an approximate manner using finite differences to represent the derivatives of Eqs. (7) and (8) on a coarser grid. The resistivities were calculated as the average value between those corresponding to bulk and interface compositions and temperatures. The nonlinear system of algebraic equations that represents the discretization of (18)-(20) was solved using MATLAB's fsolve. This was the procedure used to solve the complete distillation column model in Section 5.1. The flowchart for solving the BVP is shown in Fig. 3.

Consistency checks
In order to check the consistency of the model and the implementation, the entropy production was calculated from Eq. (35) and compared to the values obtained by Eqs. (36) and (37) for different compositions, temperatures, film thicknesses (in the range of 0.01-1 mm) and considering the effects of heat and mass coupling and interface transfer coefficients. The maximum relative error was on the order of 10 À6 , which was similar to the inaccuracies in the calculation of the partial molar entropies. This strongly suggests that the model is thermodynamically consistent. Furthermore, a careful inspection of the energy and mass balances for different parts of the column, together with a thermodynamic consistency at the interface level, also served as a check of a consistent and correct implementation of the sub-models.

A comment of the solution of the film-interface-film system in previous works
The procedure used to solve the film-interface-film system here differs from that used in the studies by van der Ham et al.
(2010) and Mendoza and Kjelstrup (2011). In previous works, they searched for the liquid film's thickness that gave the same entropy production as Eqs. (35) and (36). We argue that the representation should be valid for any choice of film thickness, which is the case for the model presented here.
In the work of van der Ham et al. (2010), the solution of the equations at the interface was obtained by discretization of the differential equations in finite volumes and transforming the resulting system using only the measurable heat flux in the vapor and evaluating all the chemical potential differences at the same temperature, i.e. the temperature of the liquid phase. In order to evaluate the differences, they assumed constant partial molar enthalpies, which introduces a small error and a thermodynamic inconsistency with respect to their choice of equation of state. The system of equations was then solved using a fixed-point iterative method, obtaining the global resistivities and the fluxes and then solving back for the conditions at each control volume. The last step implies finding a solution for a 3 Â 3 set of linear equations. This set is linearly dependent, but may be rendered nonsingular due to the numerical approximations. van der Ham (2011) reported that, for certain boundary conditions, it was impossible to find the film's thickness that satisfied the consistency condition. We believe that some of these approximations were the cause for having to search for the liquid film's thickness in order to compensate for the numerical errors.

Case studies
In order to discuss the relative importance of different transport mechanisms that are typically neglected in state-of-the-art model- ing of distillation columns for separation of oxygen and nitrogen, we shall define a set of cases. The base case (Case 0), to which the other cases will be compared, is a packed distillation column fed with a binary mixture at a rate of 10.00 mol Á s À1 (y 1 ¼ 0:79) at a temperature of 85 K and a pressure of 140 kPa. The column has structured packing Montz B1-500 (with specific packing area, a of 500 m 2 Á m À3 ), as in the experimental column in van der Ham's work (van der Ham, 2011) and its diameter was estimated at 450 mm, according to the methods in Seader et al. (2011). The column has the same packing height (2300 mm) for the rectification and stripping section, with a total area per section of 183 m 2 . The model used in Case 0 is very similar to the well-established model for distillation columns proposed by Taylor and Krishna (1992). The thicknesses of the vapor and liquid films were 0:5 and 0:1 mm respectively, which are in the typical range given by Taylor and Krishna (1992) (0.1-1 mm for the vapor film, 0.01-0.1 mm for the liquid film). The reflux ratio was set to 2 and the distillate molar flow, to 7.95 mol Á s À1 . With these conditions, the mole fraction of nitrogen in the distillate was found to be 0.987, and 95% of the oxygen was recovered in the bottoms with a purity, x B 2 , of 0.974. Five cases have been constructed to assess the effect of coupling in the films and interfacial resistivities. The details of these cases have been summarized in Table 1. All cases had the same feed conditions, operating pressure, reflux ratio, packing height, liquid and vapor films' thicknesses and distillate molar flow. Since kinetic theory provides an estimate of the interfacial resistivities with unknown uncertainty, we performed a sensitivity analysis where we multiplied all the interface transfer coefficients with the same proportionality factor, k I . This enables a broader view of the possible relevance of interfacial resistivities as they are allowed to span one to two orders of magnitude, depending on the case. In the cases without coupling, the coefficients r a iq are set to zero in the respective bulk phases.

Results and discussion
In the following, we will assess the importance of including the interfacial resistance and heat-mass coupling in the films on the overall performance of the distillation column (Section 5.1) and on the local heat and mass transfer characteristics inside the column (Section 5.2). Finally, some perspectives on the practical consequences of the findings for the modeling of distillation columns and other processes will be given (Section 5.3).  Table 1, while Fig. 5 zooms in at the stripping section. As expected, the temperature of the vapor is higher than that of the liquid, with a maximum difference of 1.5 K, where both temperatures decrease as the location moves closer to the top of the column, from about 92 K at the bottom of the column to about 80 K at the top. There is a discontinuity at z=' ¼ 0, due to the introduction of the feed. The solution of Taylor and Krishna's model (Case 0 in Table 1) and the solution of the case where coupling in the films and interfacial resistivities are incorporated into the model (Case CI) are very similar. The temperature of the vapor phase is slightly higher for Case 0, while it is lower for the liquid phase. The maximum temperature difference between Case 0 and Case CI is 0.1 K for the vapor phase's temperature and 0.2 K for the temperature of the liquid phase.

The influence of coupling and interface transfer coefficients on the overall performance of the distillation column
As elaborated in Section 1, the interface transfer coefficients, estimated by kinetic gas theory, have unknown uncertainties. The coefficients of e.g. water (Wilhelmsen et al., 2016), have been shown to deviate significantly from the predictions of kinetic gas theory. Hence, the coefficients used in Case CI should only be viewed as a reference. The interface transfer coefficients are expected to have an appreciable influence on the performance of the distillation column only if they are larger than those predicted by kinetic gas theory. In Cases CI-10 and CI-100, the interface transfer coefficients have been multiplied by a factor of 10 and 100, respectively. If the interface transfer coefficients are increased by a factor of 10, the temperature differences relative to Case 0 (where they are zero) become 0.2 K and 0.3 K for the vapor and liquid phases, respectively. We observe a significant change in the temperature profiles only if the interface transfer coefficients are increased by two orders of magnitude, where the maximum temperature difference compared to Case 0 reaches 1 K for the liquid phase, which is a large relative change compared to the temperature evolution of the fluid displayed in Fig. 4. The largest differences can be found near the bottom of the column, as shown in Fig. 5. The reason for this is that, for these operating conditions, the molar flows are smaller near the bottom of the distillation column 1 , and the differences in the calculated fluxes have a greater impact, see e.g. Eqs. (3) and (4), where the derivatives of the temperature depend on the molar flows of vapor and liquid in an inversely proportional manner.
The mole fraction of nitrogen in the vapor and liquid phases inside the column are shown in Figs. 6 and 7. The mixture is richer in nitrogen at the top of the column and in oxygen at the bottom. The mole fraction profiles display discontinuities at the location where the feed is introduced. The curves of the mole fractions for Case 0 and Case CI are nearly overlapping. Similar to the temperature profiles, this shows that the influence of coupling in the films and including the resistance of the interface is of minor importance for the overall performance of an adiabatic distillation column, provided that the interface transfer coefficients from kinetic gas theory are used. However, if the interface transfer coefficients are greater, they provide an additional resistance to mass transfer that lowers the molar fluxes. This translates into a larger nitrogen mole fraction at the bottom of the distillation column. The effect of increasing the resistances by a factor of 10 yields a maximum mole fraction difference in the vapor and liquid of 5 Â 10 À3 and 4 Â 10 À3 , respectively. If they are increased by two orders of magnitude, these differences become 34 Â 10 À3 and 29 Â 10 À3 , which are not negligible compared to the values of the mole fractions at the bottom of the distillation column.
The molar fluxes at each point inside the column (Figs. 8 and 9) are nearly identical for Case 0 and Case CI, with the relative difference being less than 1%. Even when the interface transfer coefficients are increased by a factor of 10, the molar fluxes change by less than 2%. Only when the interface transfer coefficients are multiplied by 100, the difference exceeds 10%. The area below the Table 1 Case studies and some of their main features and parameters. The parameter k I is the proportionality factor that the interfacial resistivities are multiplied by in the sensitivity analysis.

Case
Coupling curves in Fig. 8 represents the amount of nitrogen transferred from the liquid to the vapor. Since most of the curves are overlapping, the difference in the nitrogen mole fraction is very small (see Fig. 6). The measurable heat flux in the vapor is shown in Fig. 10, where the differences are larger (see Eq. (5) for the definition of the measurable heat flux). However, since the heat capacity is involved in the differential equations for the temperature (Eqs. (3) and (4)), these differences are not translated into large temperature changes.
Considering that the main differences between Cases 0 and CI are in the temperature profiles, and taking into account that the coupling coefficient in the liquid phase is large and comparable in magnitude to the heats of vaporization of the components, the explanation of the similarity in the profiles is associated with the heat fluxes and how the system adapts in order to achieve the same molar fluxes. Taylor et al. (1994) studied the sensitivity of the concentration and temperature profiles with changes in the vapor and liquid heat transfer coefficients and found that, although the heat transfer coefficients were varied over a wide range, these profiles remained similar. Since coupling affects mainly the heat fluxes, these findings are highly relevant for the present study. When we compare the interfacial temperatures for Cases 0 and CI, we find that the differences are much smaller than those between the temperatures in the respective phases. Since we are dealing with a binary system, fixing the temperature determines also the equilibrium compositions at the interface, which are very similar in both cases. Since the main mass transfer resistivities are larger than the coupling coefficients, the molar fluxes consequently also become very similar. The liquid coupling coefficient favors heat transfer from the interface to the bulk of the liquid phase, which lowers the resistance to heat transfer in the liquid film. However, the overall resistance does not change much, since the resistance in the vapor film is larger and, therefore, the rate of heat transfer is within the same order of magnitude and the tempera-   Fig. 4, zooming into the stripping section for Cases 0, CI, CI-10 and CI-100. The upper four curves correspond to the vapor, while the lower four curves belong to the liquid. ture differences between the vapor and liquid are decreased, but kept within a range. If the heat transfer rate was significantly increased, then the vapor and liquid temperatures would approach each other, while the temperatures difference would be large only if the heat transfer coefficients were small. However, if in any of these cases the predicted interface temperatures were the same, the molar fluxes would also be the same provided that there are the same compositions at each side of the interface. This is also a possible explanation of why Taylor et al.'s results (Taylor et al., 1994) are insensitive to variations in the heat transfer coefficients and why we do not observe significant differences between Cases 0 and CI in the present work. On the other hand, in Cases CI-10 and CI-100, the profiles display larger deviations since the interface provides additional resistance to mass transfer.
In summary, we find that the interface transfer coefficients must be at least about two orders of magnitude larger than the predictions from kinetic gas theory in order to have a significant influence on the performance of an adiabatic distillation column for the separation of nitrogen and oxygen. However, the films' thicknesses used in this work (Taylor and Krishna, 1992) have unknown uncertainties. If the films are thinner than those considered, then interface transfer coefficients that are smaller in magnitude could have an impact on the overall performance of the column. Further work, either experimental or with nonequilibrium molecular dynamics simulations (Badam et al., 2007;Wilhelmsen et al., 2016), is needed to reveal the true magnitude of the interface transfer coefficients, as well as the films' thicknesses.

Effects of coupling and interface transfer coefficients at the film level
In Section 5.1, we showed that the temperature and concentration profiles of an adiabatic distillation column remain to a large degree unchanged after incorporating the physical phenomena studied in this work. Part of the reason for this is likely to be the  imposed boundary conditions of the overall distillation column model, which lead the temperature and mole fraction profiles to be similar for the cases studied.
To gain further insight into how the local heat and mass transfer characteristics are influenced by the physical phenomena, we use in this section the concentration and temperature profiles from Case CI and compute the fluxes with the model configurations specified in Section 4. In particular, to calculate the fluxes across the film-interface-film system, we use the values corresponding to each location inside the distillation column from Case CI, y 1 z ð Þ; x 1 z ð Þ; T V z ð Þ and T L z ð Þ as fixed boundary conditions. This provides a range of boundary conditions, where the nitrogen mole fraction ranges from y 1 ¼ 0:09 to y 1 ¼ 0:99 and the temperatures vary from 80 to 92 K. Hence, in Section 5.1 we solved the differential equations defining the complete distillation column, but in this section we use a fixed set of boundary conditions at each location z that are the same for all the models defined in Table 1. The purpose of this study is to elucidate precisely how the physical phenomena investigated in this work influence the heat and mass fluxes, provided that the boundary conditions of the film-interface-film system are fixed.
The nitrogen and oxygen molar fluxes are shown in Figs. 11 and 12. The difference between Case CI and Case 0 comes primarily from coupling in the films. This can be inferred by observing the overlapping curves of Case C and Case CI, where the very small difference comes from the interfacial resistance. The average relative difference between the molar fluxes from Case 0 (red dotted line) and Case CI (solid line) is dramatic, 45% for nitrogen and 8% for oxygen. At some locations, the relative difference is significantly larger. Moreover, at conditions near z ¼ 1 (the top of the distillation column), the signs of the nitrogen fluxes for these two cases are the opposite of one another. Whereas nitrogen is condensing in Case CI, it evaporates in Case 0.
When the heat-mass coupling in the vapor phase was set to zero (r V iq =0), the differences with Case CI were below 1%. Hence, the main difference between Case 0 and Case CI comes from the  coupling between heat and mass transfer in the liquid film. This is expected, since the heat of transfer in the liquid phase is two orders of magnitude larger than in the vapor phase, which makes the coupling coefficient in the liquid phase larger. By comparing Case CI (red dashed line) and Case CI-100 (yellow solid line) in Figs. 11 and 12, we see that the magnitude of the nitrogen flux increases, whereas the magnitude of the oxygen flux decreases as the interfacial resistance becomes more important. Coupling of heat and mass transfer in the films has the opposite effect on the molar fluxes. Hence, both the heat-mass coupling in the liquid film and the interface transfer coefficients can change the relative rates at which the two constituents are transferred across the film-inter face-film system.
Figs. 13 and 14 show the measurable heat fluxes in the vapor and liquid phases. The measurable heat fluxes determine how much heat enters the respective bulk phases during a nonequilibrium situation. For the measurable heat flux in the vapor phase, the difference between the cases is smaller than the molar fluxes in relative terms, especially for boundary conditions corresponding to z=' > 0. By comparing the curves from Cases 0 and CI in Fig. 13, we see that neglect of coupling in the films and the interfacial resistivities leads to an underestimation of the magnitude of the measurable heat flux in the vapor phase. The coupling coefficients in both phases affect the computation of J 0;V q , where neglecting r V iq leads to 3-6% deviation and neglecting r L iq leads to 3-30% deviation.
The variable that is most sensitive to the physical phenomena studied in this work is the measurable heat flux in the liquid phase, J 0;L q , shown in Fig. 14. We observe that neglect of coupling in the films and the interfacial resistances (Case 0) yields positive values for J 0;V q . However, Case CI gives a negative heat flux, meaning that the direction of the heat flux is changed. This is due to the large value of the heat-mass coupling transport coefficient in the liquid phase, which becomes a dominating contribution to the heat flux. Large interface transfer coefficients, such as in Case CI-10 and Case  CI-100, also have a significant influence on J 0;L q , primarily in the bottom part of the distillation column. The differences between the values for each case can be very large, which indicates that J 0;L q has a large uncertainty.
For the locations where z=' < 0; J 0;L q is always negative for Case CI, always positive for Case CI-100, while, for Case CI-10, it exhibits a mixed behavior. In order to gain further insight into this, we shall study in detail the location z=' ¼ À1 (labelled A), which corresponds to the conditions at the bottom of the distillation column. The boundary conditions at this location are y 1 ¼ 0:09; x 1 ¼ 0:08; T V ¼ 92:54 K and T L ¼ 91:84 K. At Location A, the temperature in the films is shown in Fig. 15, the chemical potential of oxygen is shown in Fig. 16 and the fluxes are provided in Table 2. Fig. 15 shows that the temperature at the interface is lower than the temperature in the liquid phase. This indicates that if coupling in the films is neglected, J 0;L q should be positive, i.e. be directed from the liquid-side to the vapor-side. Since the temperature difference at the liquid-side of the interface is not very large in Case CI, coupled transport dominates over the Fourier-type transport for heat transfer, leading to a negative measurable heat flux at the liquid side. As the interfacial resistivities increase in magnitude, the figure illustrates that the temperature difference at the liquid-side of the interface becomes more pronounced, as in Case CI-10 and CI-100. With these conditions, J 0;L q is dominated by the contribution from the main transport coefficient. Larger interface transfer coefficients lead to larger temperature jumps at the dividing surface and larger temperature differences at each side of the interface. Depending on the boundary conditions, these may not be large enough to compensate for the bulk coupling effects. This explains why J 0;L q is both positive and negative depending on the location for Case CI-10.
In addition to a jump in temperature, the interface transfer coefficients also introduce a chemical potential/concentration jump at the dividing surface, which has been shown for Location A in Fig. 16. For Case CI-100, this jump (which is comparable to the variation in the chemical potential across the films) leads to a large change in the molar flux of oxygen. However, it does not change the sign of the flux, as opposed to J 0;L q .  Both coupling in the films and interface transfer coefficients have a large impact on the calculation of J 0;L q , with errors exceeding 100% and even influencing the direction of this flux. Accurate representation of this quantity may be of importance e.g. in the optimization of diabatic distillation column for increased energy efficiency, i.e. distillation columns where heat is supplied at different temperatures from the exterior. Figs. 15 and 16 show that the interface plays an important role provided that the changes in temperature and chemical potentials across the interface are comparable to those across the films. If this is not the case, it is safe to assume that the resistance to heat and mass transfer posed by the interface can be neglected.

Practical consequences for the modeling of distillation columns and other processes
The previous analysis has shown that heat and mass coupling coefficients in the bulk phases and interfacial resistivities can have a large effect on the calculation of the fluxes locally. However, these physical phenomena have little influence on the local temperature and concentration profiles of an adiabatic distillation column for separation of nitrogen and oxygen with the assumptions used in the present work. There are large uncertainties associated with both the films' thicknesses and the interface transfer coefficients. The variable that is most sensitive to the physical phenomena investigated was the measurable heat flux in the liquid phase, which can be greatly over/underestimated and even change sign. This could be a source of uncertainty in optimization studies of diabatic distillation columns, i.e. distillation columns where addition/ removal of heat to/from the liquid phase is essential. For example, Kingston et al. (2020) minimized the entropy production in a diabatic column to separate nitrogen and oxygen mixtures by searching for the optimal temperature profile of the utility. They found that, to attain that minimum, heat had to be added or removed at a specific rate at each point in the column; this rate was controlled by exchanging heat between the utility and the liquid phase. When the column is operated in a diabatic way, Eq. (4) incorporates the heat flux with the utility. Therefore, Eq. (4) contains two terms, the measurable heat flux in the liquid and the heat  exchanged with the utility. Since heat and mass transfer coupling can change the direction of J 0;L q , the heat provided by the utility can vary greatly. This leads to variations or uncertainties in the total heating and cooling duties.
When we analyze the overall effect on the performance of an adiabatic distillation column, the variations are small. This supports the current practice of neglecting coupling between heat and mass transfer in films and the interfacial resistance in steady state modeling of distillation columns as described by Taylor and Krishna (1992). A caveat is that the interface transfer coefficients could potentially be much larger than those predicted by kinetic gas theory, as their uncertainty is unknown.
Since the main variables of interest of the column (for example, condenser and reboiler duties, reflux ratio, packing height, distillate flow rate and purity) change very little, our work indicates that coupling and interface resistances can safely be neglected in process design, at least for adiabatic distillation columns operating with the oxygen-nitrogen mixture. This is in agreement with the results from Taylor et al. (1994), where they analyzed the effect of heat transfer coefficients and found that the temperature and concentration profiles did not change significantly within a certain range. This shows that even if the heat transfer between the phases is not accurately predicted through the column, one may still obtain precise results for the process design. For other types of columns where heat transfer plays a more important role (for example, diabatic columns, HIDiC), this might no longer be true and a more accurate prediction of heat transfer may be needed.
Coupling between heat and mass transfer in the liquid film influence local fluxes significantly, even though the effect on the overall performance of the column is small. With the insight from the example considered in the present work, we expect this effect to be of higher relevance in dynamic operation of distillation columns or other dynamic processes where the evolution of the state variables is less constrained by the boundary conditions.

Conclusions
The resistance to heat and mass transfer posed by the vaporliquid interface and coupling between heat and mass transfer in the adjacent vapor and liquid films is usually neglected in the modeling of e.g. distillation columns. In this work, a detailed film-inter face-film model for the transport of heat and mass from vapor to liquid has been incorporated at each spatial location in a onedimensional, continuous mathematical model of a distillation column for separation of oxygen and nitrogen operating at steady state. This multi-scale model facilitates a discussion of the influence of different physical phenomena on the overall performance of the distillation column as well as on the local heat and mass transfer characteristics.
Coupling between heat and mass transfer was shown to have a strong influence on the local behavior in the distillation column, where the effect could even alter the direction of the measurable heat flux in the liquid phase. Furthermore, the effect changed the molar flux of nitrogen by 45% on average, and the relative difference was significantly larger at some locations. This is likely to be of importance in optimization studies of diabatic distillation distillation columns. Both the heat-mass coupling in the liquid film and the interface transfer coefficients were found to change the relative rates at which the different constituents were transferred across the film-interface-film system.
The interface transfer coefficients and coupling in the films were shown to have a negligible influence on the local temperature  and concentration profiles of an adiabatic distillation column for oxygen-nitrogen separation operating at steady state, unless the interface transfer coefficients are about two orders of magnitude larger than those predicted by kinetic gas theory. Since both the interface transfer coefficients and the film thicknesses have unknown uncertainties, further studies, either with experiments or non-equilibrium molecular dynamics simulations, are needed to reveal the true relevance of the interface. The preliminary findings support the current practice of neglecting these physical phenomena in the modeling of adiabatic distillation columns operating at steady state for the separation of nitrogen and oxygen binary mixtures. Based on the strong local influence of these physical phenomena, they are expected to be of higher relevance in the modeling of dynamic processes that are less constrained by the boundary conditions.
CRediT authorship contribution statement

Declaration of Competing Interest
The authors declare that they have no known competing financial interests or personal relationships that could have appeared to influence the work reported in this paper.