In the present work, the macroscopic governing equations governing the heat and mass transfer for a general multicomponent system are derived via a systematic nonequilibrium thermodynamics framework. In contrast to previous approaches, the relative (with respect to the mass average velocity) component mass fluxes (relative species momenta) and the heat flux are treated explicitly, in complete analogy with the momentum flux. The framework followed here, in addition to allowing for the description of relaxation phenomena in heat and mass transfer, establishes to the fullest the analogy between all transport processes, momentum, heat, and mass transfer, toward which R. B. Bird contributed so much with his work. The inclusion of heat flux-based momentum as an additional variable allows for the description of relaxation phenomena in heat transfer as well as of mixed (Soret and Dufour) effects, coupling heat and mass transfer. The resulting models are Galilean invariant, thereby resolving a conundrum in the field, and always respect the second law of thermodynamics, for appropriate selection of transport parameters. The general flux-based dynamic equations reduce to the traditional transport equations in the limit when mass species and heat relaxation effects are negligible and are fully consistent with the equations established from the application of kinetic theory in the limit of dilute gases. As an added benefit, for the particular example case of hyperbolic diffusion we illustrate the application of the proposed models as a method to allow the use of powerful numerical solvers normally not available for solving mass transfer models more generally.

## I. INTRODUCTION

Considerable effort has been devoted already in the macroscopic modeling of mass and heat transfer within multicomponent mixtures. Regarding multicomponent diffusion, one can trace its origin to a generalization of the empirically derived Fick's law for a binary system.^{1–3} However, for a more formal description, one had to wait the development of nonequilibrium thermodynamics (NET).^{2} This also allowed one to take advantage of the Onsager–Casimir's reciprocal relations^{4–6} that resulted in significantly reducing the number of independent parameters, *n*(*n*-1)/2 instead of (*n*-1)^{2},^{2} where *n* is the number of diffusing species. In parallel with Fick's work, Maxwell developed an inverse relation correlating the gradients of concentrations to diffusive fluxes^{7,8} that was later extended to multicomponent systems by Stefan.^{9,10} This fundamental relation was eventually extended to include pressure and thermal effects by Curtiss and Hirschfelder^{11} and was finalized shortly afterward in the monumental treatise of “molecular theory of gases and liquids” by Hirschfelder, Curtiss, and Bird.^{12} More recently, the two approaches (the generalized Fick and the generalized Maxwell–Stefan) were shown to be related in the work of Curtiss and Bird^{13}—in a recent article by Bird and Klingenberg^{14} these and other historic instances related to multicomponent diffusion have been nicely outlined and summarized. Of interest also is the very recent article by Swaney and Bird,^{15} where through the use of thermodynamic relations the authors develop equations for energy, mechanical energy, and entropy for multicomponent mixtures.

In parallel to the above works, several efforts have been dedicated to obtaining more general expressions for multicomponent diffusion applicable to arbitrary fluid systems. From the early works, it is worth noting the work of Truesdell^{16} in developing a more general form of Maxwell–Stefan theory using general principles of the mechanics of arbitrary continua. However, these efforts in modeling multicomponent diffusion within the at the time available continuum mechanics formulation led to indeterminate results when, by theoretical necessity, the modeling of the rate changes of the relative species momenta was also included.^{17} This unsatisfactory situation was specifically mentioned by Truesdell and Noll in Sec. 130 of their fundamental work, the “non-linear field theories of mechanics” and that statement remained present even in the last edition.^{17} In retrospect, as part of the outcomes from the present investigation, we attribute that problem in the introduction of species “stresses” when the balance of species momenta was introduced in that theory.^{17} We feel that species-specific stresses cannot be physically defined from a continuum point of view, as the individual species are not isolated in individual parts of space but are sharing the same space with all others. As a result, any physically imposed constraints following that definition are invalid. This problem issue is common to all such continuum theories that make use of similar assumptions, usually known as standard “mixture” theory. Later work by Lam^{18} trying to resolve the indeterminacy issue at the continuum level by resorting to collision theory-based explanations led to formulas too complicated to use (on top of the need to use several questionable assumptions). Other works, when attempting to generalize the multicomponent diffusion model, they did that at the constitutive level of fluxes, generalizing the Stefan–Maxwell relations through the use of nonequilibrium thermodynamics at that level—see, for example, Refs. 19–21.

The standard mixture theory and/or the use of species stresses has continued to be used in many recent works that have attempted to extract better expressions for the generalized Stefan–Maxwell relations.^{22–27} The difference between those works and the more traditional ones mentioned at the end of the previous paragraph that sought to express the diffusion fluxes statically in terms of constitutive relations directly relating them to potentials (such as gradients of the concentration and temperature or the chemical potential) is that in those an additional attempt is made to model their dynamics, following a multi-fluid approach. Such an approach has been followed in several other works in the literature, (primarily the two-fluid case; see, for example, Refs. 28–32) albeit in those works inertia effects and flux dynamics were ultimately neglected. Instead, the emphasis in those works was to capture the effect of nonequilibrium stresses to the flow developing in one of the two-fluid phases, typically corresponding to a polymeric system^{28–31} but also for suspensions of solid particles.^{32} Of special importance is previous work^{33–35} where a non-equilibrium thermodynamics-based, two-fluid theory was presented to describe the structure and rheology of shear banding rodlike micellar solutions. Even if inertia in diffusion fluxes may be often negligible (and a recent work, based on a statistical analysis near equilibrium, claims that this may be what happens for collision-dominated dilute ideal gases^{36}) developing a consistent theory that can account for it may be very worthwhile. First, it is needed to predict finite speed of propagation of diffusive processes within matter, a phenomenon for which there is considerable experimental evidence,^{37} thus removing the aphysical infinite speed characterizing classical descriptions of heat and mass transfer. Second, it establishes needed consistency to the continuum modeling of transport processes and third it may be the best way to bring in the connection to nonequilibrium thermodynamics—it may also have significant applications under conditions far away from equilibrium as well as it may lead to novel numerical approximations even when flux inertia effects can be safely neglected (an example of which is illustrated further in this work).

In the present work, we elaborate on the modeling of flux dynamics for simultaneous heat and mass transfer in multicomponent mixtures with the objectives (a) correcting the mixture theory by consistently including inertial (kinetic energy) effects to diffusion and heat fluxes that may be important under conditions where an appreciable differential motion between various populations of the chemical species is present, (b) developing the theory for diffusion and heat transfer at the same level (using suitably defined momentum densities) as for momentum transfer. This is worthwhile as it offers a much needed closure to the full equivalence between all transport phenomena to which R. B. Bird has so much contributed, and (c) demonstrating a more straightforward evaluation of the extraction of generalized Stefan–Maxwell relations, when inertia flux effects can be neglected, as their connection to nonequilibrium thermodynamics becomes more straightforward at the momentum densities level. The addition of a momentum density associated with the heat flux allows one to consistently describe not only relaxation effects in heat transfer^{38–40} but also the mass and heat transfer coupling (the Soret and Dufour effects).

To achieve these objectives, we employ the single generator bracket formulation of non-equilibrium thermodynamics (SGBF-NET).^{41} This represents one of the simplest formulations of the most modern implementations of nonequilibrium thermodynamics, the GENERIC formulation^{42} to which it has been shown to be completely equivalent when applied to macroscopic modeling.^{43} The applicability of SGBF-NET and GENERIC to model transport problems has been shown many times, both in a variety of standard transport cases^{41,42} but also in many special cases, most recently that corresponding to developing thermodynamically admissible 13-moment equations, as described in another article belonging to the present special *Physics of Fluids* issue.^{44} Moreover, it is a thermodynamically consistent formulation that as we show can reproduce the exact functional forms for the heat and mass fluxes obtained from the kinetic theory of dilute gases.^{12} In this way, for this particular example, not only one can have full expressions (that are non-trivial!) for the dependence of the phenomenological transport coefficients on the system variables, but also one can hope, through the use of minimal adjustable parameters that can in principle be determined experimentally, to be able to extend these developments to more complex cases (e.g., dense gases and/or liquids) where any *a priori* kinetic theory-based prediction is lacking. However, we also need to note here that one may have a variety of ways through which heat transfer is carried over, depending on the material of interest and the conditions. The form used in the present work is suitable to describe heat transfer as established through collisions and based on the kinetic energy of the participating molecules, exactly as assumed in the kinetic theory of gases.^{12} Other forms of heat transfer may be possible, necessitating, in general, separate internal variables and fluxes—see, for example, Refs. 38–40 and 45.

In Sec. II, we show how the SGBF-NET can be used to naturally extend the modeling of heat and mass transfer within multicomponent systems to also include the effect of possible heat and different species momenta. In Sec. III, we show how this treatment is compatible in the limit of small differential velocities, to classical descriptions of heat and mass transfer, including the Fourier heat conduction, the Fickian diffusion in a binary mixture, and the Maxwell–Stefan equations, also corrected for coupling with heat transfer (Soret/Dufour effects). In fact, related to the latter case, we also show in Sec. III how the present development leads to equations that are fully compatible to those obtained from kinetic theory for dilute gases^{12} but also allowing to consistently extent those equations to other systems. This analogy allows us in the dilute gases case to obtain specific values for all the transport coefficients entering the theory. This is demonstrated, for a specific binary gas case, also shown in Sec. III. One further specific example is examined in detail in Sec. IV: That of the binary diffusion of a dilute solute in a liquid solvent. We show there, beyond the compatibility of the flux model to the standard Fickian diffusion in the limit of small differential velocities, the usefulness of the flux approach in setting up very simple but also very powerful numerical approximations, exemplified in the simple problem of transient diffusion in a slab. More specifically, taking advantage of the hyperbolic character of the flux-based model equations, we show how the use of explicit time integration solvers can reproduce results that otherwise would have required significantly more computational effort, demanding implicit numerical time integrations. Finally, in Sec. V, we present our conclusions.

## II. SGBF-NET THEORY AND DERIVATION OF MULTICOMPONENT TRANSPORT EQUATIONS

We apply here the single generator bracket formalism of nonequilibrium thermodynamics (SGBF-NET) as developed in Ref. 41. As shown in Ref. 43, this is equivalent for macroscopic systems to the more general GENERIC approach.^{42} For completion, we start with a very short presentation of the bracket formalism of nonequilibrium Thermodynamics, borrowing from the description offered in a recent publication^{46} that outlines its influence in elastic fluid modeling—for a more detailed presentation the interested reader is referred to the above-mentioned research monograph.^{41}

The bracket formalism offers a generalization of the Poisson bracket description (which is the most general description of Newtonian mechanics, see, for example, Refs. 47 and 48) from conservative, i.e., Hamiltonian, to dissipative systems. Of critical importance was the development first of Poisson brackets for conservative continuum systems, such as the Euler equations, by Glebsch^{49} and Arnold.^{50} The development of the dissipative bracket followed next.^{51–53} The bracket formalism also extends equilibrium to non-equilibrium thermodynamics, with key assumption that there is a set of variables (superset of the equilibrium ones) in terms of which a local equilibrium can be defined, one of the variables being a local entropy density. Within the bracket formalism, there are four steps that lead to the development of the system of evolution equations governing any given system:^{41,46}

First, one needs to define a set of acceptable state variables, ${z1z2\u2026zn}$, in terms of which the system is specified. Those are typically field densities, including, at a minimum, those needed to specify the system at equilibrium.

Second, the Hamiltonian needs to be specified, physically identified with an (extended) energy of the system, as a functional (integral over the system volume) of the state variables, $H=\u222bVh(z1z2\u2026zn;\u2009\u2207zi;\u22072zi;etc.)dV$, the dependence on the first, second, and higher gradient being optional.

Third, the Poisson bracket (governing the reversible component of the system dynamics), ${F,G}$, needs to be specified. It has all the properties of an acceptable Poisson bracket used to specify conservative Hamiltonian dynamics, i.e., it has to be an antisymmetric bilinear functional of the Volterra derivatives for any pair of arbitrary functionals,

*F*,*G*, that also satisfies the Jacobi identity.^{41}Fourth, the dissipation bracket (governing the irreversible component of the system dynamics), $[F,G]$, needs to be specified. This last component is introduced to account for dissipative phenomena; it is linear in terms of the first functional, but in general nonlinear in terms of the second. It also needs to satisfy (upon linearization, close to equilibrium) certain positive definite symmetric or antisymmetric properties between the various Volterra derivatives of the functionals of

*F*,*G*depending on the symmetry upon time inversion of the variables with respect to which the Volterra derivatives are taken.^{41}These properties guarantee satisfaction of the Onsager/Casimir reciprocal relations of linear irreversible thermodynamics^{4–6,41}as well as the satisfaction of the second law of thermodynamics.

Once those four steps are carried out, the dynamics of any arbitrary functional *F* are simply specified as

If so wished, the evolution equations for the Eulerian partial time derivatives of anyone of the state variables, $\u2202zi\u2202t$, can then easily be deduced by comparing Eq. (1) against the expression obtained through the application of differentiation by parts

and equating the factors in front of each Volterra derivative $\delta F\delta zi$.^{41} The key advantage of the bracket approach lies in (a) being able to separate cleanly reversible and irreversible (dissipative) components in the dynamics by simply tracing their provenance to the Poisson and dissipative brackets, respectively, and (b) automatic satisfaction at the macroscopic level of all the irreversible thermodynamic requirements.

In the following, we will apply the bracket formalism to the flux-based modeling of heat and mass transfer. The first step involves the selection of the appropriate state variables. To describe heat and mass transfer between *n* species, $aj$, where $j=1,\u2026,n$, the standard set of variables are: (a) each one of the *n* species mass densities, $\rho j$, and (b) their corresponding momentum densities

where $v\xaf(j)$ is the jth species mass-averaged velocity, (c) the entropy density, *s*. Note that the total mass density, *ρ*, and its corresponding total momentum density

(where the second equality is used for the definition of mass averaged velocity, $v\xaf$) are then obtained through a straightforward summation of the corresponding species quantities

For reasons that will become obvious in Sec. III, we propose to consider here the following equivalent set of primitive variables: (a) the species mass densities $\rho j$ [through which the total mass density *ρ*, can also be obtained, from the first relationship in Eq. (5)] (b) the total momentum density, $u\xaf$, (c), the reduced momentum densities corresponding to each species

(d) the total entropy density, $s$. Note that in this case, the total momentum density ($u\xaf$) becomes an independent primary variable of the system, reestablishing the status it enjoyed for single component system. However, in exchange, the reduced momentum densities are not all independent. Instead, they must satisfy the following constraint:

In addition, note that, to ensure symmetry and consistency between the treatment of heat and mass transfer, we introduce one more variable, the reduced heat momentum density $u\u2322\xaf(0)$. This is defined in terms of the heat flux $q\xaf$ as

where the effective heat mass density, $\rho 0$, is a material property with units of mass density entering in the definition of the contribution to the extended free energy density due to the heat flux (see below), and $v\u2322\xaf(0)$ is a corresponding to the heat flux velocity with the physical meaning that it would have resulted to the same overall increase in the entropy if it was added to the fluid velocity. It therefore plays the same role for heat transfer as the relative species velocity $v\u2322\xaf(i)$ does with respect to the diffusive flux of species *i*. After the addition of this new reduced heat momentum density, the appropriate state $2n+3$ variable space can be defined as $w\u2322\xaf\u2261[s,\rho 1,\u2009\u2026\u2009,\rho n,u\xaf,u\u2322\xaf(0),u\u2322\xaf(1),\u2026\u2009,u\u2322\xaf(n)]T$, with the understanding that not all the reduced species momenta are independent, as they need to obey the homogeneous constraint described by Eq. (7), thus effectively reducing the dimensions of the state space to $2n+2$.

The second step in the bracket formulation involved the description of the total energy (Hamiltonian) of the system. This is provided here as

where the last term is written in terms of the prescribed system variables with $\epsilon =\epsilon (s,\rho 1,\u2026,\rho n)$ representing the equilibrium contribution to the energy density. Note that the heat mass density is not an independent state variable, but rather, an additional material property, $\rho 0=\rho 0(s,\rho 1,\u2026,\rho n)$. At this point, we cannot be more specific as to what this effective mass density is; from a units point of view, one could alternatively express it in terms of the entropy density, s, the temperature, T, and a characteristic speed, $cs$ (in an equation, that is, introduced for convenience and shall be used as a proper definition for the characteristic speed, $cs$) as

The additional term correcting the total Helmholtz free energy (at the lowest order and taken here as isotropic) for the presence of nonzero heat fluxes appears in Eq. (9) similar to the term introduced by Coleman *et al.*^{54} with the coefficient $z=q0(sT)2$ representing the magnitude of an isotropic tensor $Z\xaf\xaf=zI\xaf\xaf$ in the corresponding general tensorial and anisotropic expression proposed by Coleman *et al.* [Eq. (1.47) in Ref. 54].

Based on the general expression for the Hamiltonian provided by Eq. (9), one can readily calculate the corresponding Volterra derivatives of the Hamiltonian as

where, as usual, the jth chemical potential and temperature are defined as $\mu j\u2261Mj\u2202\epsilon \u2202\rho j$ and *j* takes the values $j=1,2,\u2026,n$. Note that the Volterra derivatives with respect to the reduced species momentum densities had to be corrected in order to satisfy the same constraint as the reduced momentum densities, Eq. (7) (which can be expressed as the *n*-vector generated by the reduced momentum densities to be normal to the $[11\cdots 1]T$ (*n* elements) row vector). Note that there is one-to-one correspondence between the set of variables used here and that involving the species momenta, with their relationships as well as that between the corresponding Volterra derivatives of an arbitrary functional *G*, summarized for convenience in Appendix A. *The important feature of the variables chosen as system variables above is that now the Volterra derivatives with respect to the species densities have their unique contributions arising from the species kinetic energies that are Galilean invariant*, i.e., *independent of the velocity of the coordinate system in terms of which we represent our system.*

Now, if desired, we can also involve additional state variables, for example, as in our previous work,^{33–35} we can associate with each component *k*, also a conformation tensor $C\alpha \beta (k)=\u27e8\rho kR\alpha (k)R\beta (k)\u27e9\u2261\rho kc\alpha \beta (k)$ as the mass-weighted average of the second moment $c\alpha \beta (k)$ of a structural position vector, $R\alpha (k)$, corresponding to component *k* (assumed to be a macromolecule). Note that in such a case, there is an extra contribution to the Hamiltonian through an “elastic” free energy density, $he=he(\rho k,C\xaf\xaf(k),T)$ that also contributes to the Volterra derivative $\delta H\delta \rho k$ an extra elastic energy related term $\u2202he\u2202\rho k$. This effect can lead to very interesting phenomena, such as stress induced migration,^{28} and as such will be studied in a latter work; however, for simplicity it is not considered in the following work presented in this manuscript.

Third, the Poisson bracket, which describes the only reversible process (convection) taking place, can only involve one and only one genuine velocity field, $v\xaf$, corresponding to the total momentum, $u\xaf$. This is in contrast to the multi-fluid model.^{35} The reduced momenta, $u\u2322\xaf(i)$ must now be considered as internal variables, and as they are vector quantities need to obey either a covariant or contravariant dynamics to lead to materially objective equations. In this case, the covariant formalism is the appropriate one (as the fluxes are ultimately driven from the gradient of scalar potentials) leading to a Poisson bracket that is cast as follows:^{41,46}

where the Einstein's summation convention over repeated Greek indices (from 1 to 3) is followed here and in the following. The covariant choice is also confirmed later when comparing the results to those of kinetic theory.

Fourth, as the processes of heat transfer and mass diffusion, which are an inherently irreversible process, are driven by the reduced momenta, it is in the dissipation bracket where new terms involving these new variables will appear, replacing previous contributions described directly in terms of the gradients of temperature and the chemical potentials. Therefore, this is where the extension to the bracket formulation of the standard evolution of any arbitrary functional $F$ in a multicomponent system, such as for example, described in Sec. 7.3 of Ref. 41, is necessary. The specification of the dissipation is facilitated upon introducing its primitive component. The primitive dissipation bracket, $[F,H]p$ (defined as the dissipation bracket without the entropy correction term), can then be inferred as consisting of two terms, one antisymmetric (therefore leading to no dissipation and entropy correction), $[F,H]a$, and a symmetric one, $[F,H]sp$

The antisymmetric dissipation bracket is cast so that it faithfully restores the individual species mass density evolution equations to their standard form when the reduced momenta are described in terms of the species mass average velocities [as duly confirmed below—see Eq. (22)]. Similarly, it also introduces the convective effects of the diffusive fluxes to the entropy flux as well as the entropy flux due to the heat flux. The contributions of the species reduced momenta coupling with the species mass balance equations are formally generated from the Poisson bracket that would have been created if we could have multiple species velocities [first term on the right-hand-side of Eq. (12)] with the mass average velocity $u\xaf$ being replaced by the species velocities $u\xaf(i)$ and then the Volterra derivatives with respect to the species velocities being expressed in terms of the appropriate elements of the variable space $w\u2322\xaf$, i.e., the reduced species momenta $u\u2322\xaf(i)$ and the mass average velocity $u\xaf$ [see Appendix A, Eq. (A3)].

The contributions of the Volterra derivatives with respect to the mass average velocity $u\xaf$ are the only ones that can be accommodated as part of the Poisson bracket, Eq. (12), based on the convective, reversible dynamics, structure that they represent. We therefore need to allocate the remaining contributions to the dissipation bracket, $[F,H]a$, representing what constitutes the component of the antisymmetric portion of the dissipation bracket, $[F,H]a$, that is, related to the coupling of the species reduced momenta with the mass densities. In addition, as the diffusive fluxes also transfer entropy, additional terms need to be added coupling the species reduced momenta to the entropy density. Finally, an equivalent term for the heat flux reduced momentum needs also to appear to model heat transfer. When all these contributions are included, we find

The primitive part of the second, symmetric, contribution of the dissipation bracket, $[F,H]sp$, follows from the general form, which is assumed to be a general bilinear functional with respect to *F, H,* near equilibrium, as follows:

where $Q\alpha \beta \gamma \epsilon \u2261\mu (\delta \alpha \gamma \delta \beta \epsilon +\delta \alpha \epsilon \delta \beta \gamma )+\kappa \u2032\delta \alpha \beta \delta \gamma \epsilon $ with $\mu ,\u2009\kappa \u2032$ are the shear and bulk viscosity^{41} and $\alpha ij,\u2009\alpha i0,\alpha 00$ are the *i-j* generalized binary mass transfer coefficients, $i,j=1,2,\u2026n,j\u2260i$, *i*th generalized Soret/Dufour coefficients, $i=1,2,\u2026,n$, and a scalar generalized heat conductivity coefficient corresponding to this problem, respectively.

As shown in Appendix B, the expression offered in Eq. (15) is the most general which (a) is bilinear in terms of the Volterra derivatives involved and (b) only considers an isotropic behavior (thus involving only scalar transport coefficients). The latter assumption, reasonable in the case that there are no other, tensorial in nature, structural variables, may not hold when those are present, as for example in the case of polymeric systems—see, for example, Ref. 55. For simplicity, we keep it here with the caveat just noted. Moreover, please note that by construction only the symmetric, off diagonal part of the matrix of the coefficients $\alpha ij$ matters, and these are the only ones considered [for a total of *n*(*n*-1)/2 independent components]. Furthermore, note that not all the generalized Soret/Dufour coefficients are independent as, due to the fact that the reduced species momenta are required to satisfy the homogeneous constraint defined by Eq. (7), they can be defined only up to a constant. To leave them fully specified, the Soret/Dufour coefficients have therefore to be required to satisfy the corresponding homogeneous (conservation of linear momentum) constraint

In addition, it is important to note that as a result of the symmetry of the transport coefficients $\alpha ij$ and the constraints described by the above equation, the Onsager/Casimir reciprocal relations are duly satisfied. Moreover, for the non-negative requirement for the corresponding entropy production, all, that is, required to ensure that the quadratic expression for the local rate of entropy production (see below for the derivation) is positive semidefinite. This generates straightforward but rather complex general requirements on the transport coefficients $\alpha ij,\u2009\alpha i0,\alpha 00$; suffice to mention though here that in the particular case of a binary mixture those requirements reduce to (a) the transport coefficients $\alpha 12,\u2009\alpha 00$ to be non-negative and (b) the Soret/Dufour coefficients $\alpha 10=\u2212\alpha 20$ to satisfy the inequality $(\alpha 10)2\u22642\alpha 12\alpha 00$. We also note that general constraints can be obtained through the use, for example, of a modified Sylvester's criterion: all principal minors of the corresponding square matrix need to have non-negative determinants.^{56}

The entropy correction to the above dissipation bracket is calculated, as usual, based on the constraint of the conservation of the overall system energy (Hamiltonian) *H*, as follows.^{41} In general, if the primitive dissipative bracket is offered as

then, the entropy correction term, $[F,H]ec$ is defined as

with the (full) dissipation bracket, $[F,H]$, provided by the sum of the primitive and entropy correction terms

Corresponding to the above-described Poisson and dissipative brackets, the evolution equations for the state variables are then developed as usual,^{41} i.e., by requiring the total derivative of any arbitrary functional $F$, also expressed as [see also Eq. (2)]

To be identical to the fundamental equation in the Hamiltonian bracket nonequilibrium thermodynamics, Eq. (1). Using this rationale, the expression for the total time derivative of individual component densities, $d\rho idt$, would be equal to the coefficient of $\delta F\delta \rho i$ in the one generator expression for $dFdt$, based on the evaluation of the Volterra derivatives of the Hamiltonian provided by Eq. (11), leading to

which duly reduces the equation to the familiar form of the continuity equation for a single component

Given that the mass densities are independent and additive, performing a summation over all *n* components in the system of the above equation gives us the overall continuity equation of the system

where $\rho =\u2211i=1n\rho i$.

Through similar reasoning, one can arrive at an expression for the total time derivative of the total momentum density $u\xaf$; that comes out to be

where $\sigma \beta \alpha (c)$ is a reversible contribution to the stress, as it arises from the reduced species momenta

and $p$ is the non-equilibrium pressure of the system defined as

Note that the above expression also accounts for purely nonequilibrium effects arising due to differences in the velocities of the components. However, it is independent of the overall velocity of the system (and therefore Galilean invariant)! Indeed, if one focuses on the contributions that the overall kinetic energy may have on the definition provided by Eq. (26) one can easily find those to be zero

This is an important finding and necessary for the well-posedness of the proposed theory. From Eq. (23) and the definition of total momentum, Eq. (4), one can get the following evolution equation for the mass average velocity:

Similarly, repeating the exercise for the system entropy density, *s*, we get

where $\sigma s$ represents the total rate of entropy production and is given as

As mentioned, requiring this expression to give us always non-negative results for the allowed (through constraints) variables, gives as constraints for the various transport coefficients that are involved.

Repeating the exercise for the reduced momenta densities results in the following expression:

Equation (31) can be also written in a slightly more compact form, by combining terms in the definition of the time derivative as well as those involving the Volterra derivatives with respect to the species mass densities, as

with the overscript triangle representing the generalization of a lower convected time derivative for a covariant vector quantity^{57} weighed by a density. This is a particular case of the general Oldroyd time derivative,^{46,58} defined as

Equation (32) can be further rewritten by using Eq. (11) to substitute for the Volterra derivatives of the Hamiltonian as

Note that at equilibrium (i.e., in the absence of any flow and inhomogeneities) the values of the reduced momenta densities $u\u2322i\xaf$, and their lower convected time derivatives, are equal to $0\xaf$.

Regarding the species momenta, their evolution equation is obtained from Eqs. (22) and (28) when the definition of the reduced momenta, Eq. (6), is substituted into Eq. (31)

The last equality was obtained by using the identities

Note that, as far as the explicit dependencies on the velocities are concerned, Eq. (35) duly agrees to the results from kinetic theory.^{59,60} This agreement would not have been possible if another form (i.e., such as that associated with a contravariant vector, or scalar) was used to represent the reversible (i.e., dictated by the Poisson bracket) time derivative rather than the covariant one adopted here. Moreover, in comparison of Eq. (35) with Eq. (130.4) in the last of Truesdell's voluminous work,^{17} the only difference is that we include a term $\rho i\rho (\u2207\beta \sigma \beta \alpha (c)\u2212\u2207\alpha p+\u2207\beta (Q\alpha \beta \gamma \epsilon \u2207\gamma v\epsilon ))$ in lieu of the divergence of a partial stress term $\u2207\beta T\beta \alpha (i)$ which though is left undefined except for the requirement that summing all the partial stresses results in the total stress of the system. In fact, Truesdell and Noll in the last few pages of their work, following Eq. (130.4),^{17} attempt to make sense of that term based on the work of Adkins and coworkers, with the end result making them feel unsatisfied as it is described in the last paragraph of their work:^{17} “The foregoing summary makes it clear that the rational theory of diffusion is in its infancy, with many possibilities awaiting search. The uncertainty as to even the overall pattern for the right constitutive equations, reflected by the diversities among the five papers on the subject issued by the same author within a few months, suggest that some basic principle of the subject remains to be discovered.” Thus, in addition to offering a meaningful correction, *a key benefit resulting from our investigation is an explicit evaluation of all stress, pressure, and force terms.*

This concludes the formal description of the general flux-based heat and mass transfer equations. In Secs. III–IV, we will explore different limits of these equations to (a) validate them, (b) interpret important effects included within them as subcases, and (c) explore the symmetries and mathematical characteristics in developing efficient numerical schemes.

## III. LIMITING CASE OF SMALL FLUXES: VALIDATION AND GENERALIZATION OF STEFAN–MAXWELL AND SORET/DUFOUR EQUATIONS AND AN APPLICATION TO A BINARY SYSTEM

First, it is of interest to see what the theory predicts for the limit of small differential velocities (therefore neglecting square velocity terms) and small rates of change (therefore equating to zero the corresponding Oldroyd time derivatives). In that limit, Eq. (34) becomes

where, in obtaining the second equality in the first line, the constraint given by Eq. (16) has been used. Note that, as a result of the constraint provided by Eq. (16), in the above equation one can substitute the absolute species velocities, $v\xaf(i)$ by the relative ones $v\xaf\u2322(i)$ through which one can easily evaluate the corresponding fluxes, through Eq. (6).

These equations bear high resemblance to the Maxwell–Stefan equations, with corrections for the coupling with heat transfer (Soret/Dufour effects). The best way to check them (as well as to get some idea on possible expressions for the transport parameters entering into these equations) is through a comparison against the results of an independently validated microscopic theory, namely the kinetic theory of gases. Indeed, the original Ref. 12 provides the following two relations for multicomponent heat and mass transfer. First, Eq. (7.6-12) of Ref. 12 provides an expression for the heat flux (after correcting for the enthalpy flux included in (7.6–12) as it represents in the work the “total energy flux”), which after rewriting it slightly to our notation is given as

where $\lambda $ is the thermal conductivity, and the $Dij\u2009and\u2009DiT$ are the binary diffusion and thermal diffusion coefficients, respectively.

A direct comparison of Eq. (38) against the first line in Eq. (37) shows that they are identical when the following expressions are used for the transport coefficients $\alpha i0,\u2009i=0,1,\u2026,n$:

A physical interpretation shows the specific form of the more general transport coefficients, defined in Eq. (15), as appropriate for ideal gases. The first relation shows that the generalized heat conductivity coefficient is proportional to the entropy density squared, temperature and inversely related to the thermal conductivity. As shown in Eq. (30), it also weighs the square of the heat flux velocity in the expression for the rate of entropy production. The second set of generalized Soret/Dufour coefficients, connecting the reduced momentum densities [first term of Eq. (34) to the heat flux, they involve both the thermal diffusivities and the binary mass diffusivities]. Note that these satisfy the constraint provided by Eq. (16) and the Onsager/Casimir reciprocal relations. Importantly, mapping the general result to the specific case of an ideal gas provides explicit limiting forms for the unspecified generalized transport coefficients.

Second, Eq. (8.1-3) of Ref. 12 provides an expression for the mass transfer which after some rewriting to our notation reads as

where the vector $d\xaf(i)$ is defined in Eq. (8.1-2) of Ref. 12 that in the absence of external forces reads as

Now, for an ideal multicomponent gas, the chemical potential of component *j*, $\mu j$, is given as [see, for example, Eq. (6.1.8) of Ref. 61]

Therefore, with the use of Eq. (42), the gradient of the chemical potential, $\u2207\xaf\mu j$ at constant temperature, with temperature (or, equivalently, entropy) being an independent variable of the problem, is given as

leading to the following equality:

Using Eq. (44), and the known value for the parameter $\alpha i0$ as provided by Eq. (39), Eq. (40) becomes

If we then use the first line in Eq. (37) that we have shown above to be equivalent to Eq. (38) derived from the kinetic theory [Eq. (7.6-12) from Ref. 12] to eliminate the dependence on the temperature gradient in Eq. (45) it yields

or, using Eq. (39) to eliminate $\lambda $ and after multiplying both sides by $cRGT$ we find

Finally, this can be written explicitly in terms of velocity differences by using the condition given by Eq. (16) as

Comparison of Eq. (48) with the second line of Eq. (37) shows that the two equations are fully equivalent provided that the following expression is used for the phenomenological coefficients $\alpha ij$ for dilute gases:^{12}

Nicely, the expression is symmetric as required.

Thus, in this limit the mass transfer coefficients, which by Eq. (34) relate the reduced momentum density to the species mass velocities, explicitly depend not only on the mass diffusivities but also on the Soret/Dufour and generalized heat conductivity transport coefficients. This result is to be expected as the relative mass velocities are coupled to the heat flux. This comparison shows (a) the complete equivalence between our model and the one arising from the kinetic theory of dilute gases^{12} and (b) identifies through the determination of the entropy production (and therefore of the non-negative definite constraints on the generalized mass transfer coefficients) the useful range where this model is thermodynamically admissible.

It remains to see now some typical values for those transport coefficients. Let us first consider the subcase of an ideal binary gas mixture. First, notice there that, because of the constraint expressed by Eq. (16), $\alpha 10=\u2212\alpha 20$ (as mentioned before), and the only nonzero mass flux coefficient following Eq. (49) is $\alpha 12$, which is always non-negative:

where

For dilute binary gas, the significance of the reduced coefficient $\alpha \u0302ij$ emerges when one considers the only other criterion that we can develop for the entropy production as described by Eq. (30) to be non-negative besides the coefficients $\alpha 12,\u2009\alpha 00,\u2009\mu ,\u2009\kappa \u2032$ being non-negative

which is of course duly satisfied given the correspondence provided by Eq. (51).

It is of interest to evaluate a typical value of the correction suggested by Eq. (50) on the reduced mass transfer coefficient $\alpha \u030212$. To do so, let us use one of the many systems discussed in Ref. 12 for which a complete set of data on the binary diffusion, thermal diffusion, and thermal conductivity are available. A nitrogen–hydrogen, $N2\u2212H2$, at 65% of hydrogen binary system seems to fit that constraint: In Table 8.4–11 of Ref. 12 it is reported a thermal conductivity at T = 273.16 as $\lambda =2.0\xd710\u22124\u2009cal/(cm\u2009s\u2009\xb0K)$, whereas in Table 8.4-12 a binary diffusivity at the same temperature is reported as $D12=0.67\u2009cm2\u2009s\u22121$. Finally, from Table 8.4-14 for the same system and at about the same temperature a value for the thermal diffusion ratio can be estimated as about $kT\u22480.068$. From the definition of the thermal diffusion ratio for a binary gas mixture [Eq. (8.1-13) of Ref. 12]

Comparison of Eq. (53) with Eq. (40) applied to the same binary gas mixture gives us the following connection between the thermal diffusion ratio and the Soret/Dufour thermal diffusion coefficients:

This results in the following relationship for the $\alpha 10$ coefficient:

Finally, the total molar concentration for an ideal gas at the same temperature is $1/22\u2009400\u2009mol\u2009cm\u22123$. With all this information, it is straightforward to evaluate the relative correction to the reduced mass transfer coefficient $\alpha \u030212$

We therefore see that at least in this case this correction is small and can be neglected, thereby approximating $\alpha 12\u2248\alpha \u030212$, a practice followed in the literature and also below (whether though this is true in all cases it remains to be seen).

Another case where the equations can be tested for agreement with previous models is in the Fickian binary diffusion of a dilute solute in a liquid solvent in the absence of heat transfer and, as mentioned above, neglecting the correction to the reduced mass transfer coefficient, and using the expression appropriate for the ideal gas case. Under those circumstances, the approximations often made are

leading to the following expression for the right-hand-side and the first term of the left-hand-side of Eq. (47), respectively, as applied to that liquid binary system, for $i=A,\u2009j=B$:

where $j\xafA$ is the mass-based diffusive flux for component $A$, corresponding upon equating the two sides to the following expression for the flux:

which is the appropriate expression for Fick's law for this system [see formula (A) in Table 1.3 in Ref. 62]. It is instructive to note that under those conditions Eq. (32) becomes

A final exploitation of the results obtained here is in delivering an approximation where the fluxes are eliminated as independent variables of the problem but are rather reconstructed as proportional to the gradients of the corresponding potentials (typically followed in heat and mass transfer applications). The Hamiltonian bracket formalism allows us easily to obtain the corresponding equations following an approach that was first introduced systematically in Refs. 63 and 64 in connection to a dilute emulsion and taking as an example that of a damped linear oscillator. It simply suffices to (a) remove those variables from the Poisson expression, Eq. (12), (second and third line in that equation), (b) remove the asymmetric contributions, Eq. (14), from the dissipation bracket, and (c) write the remaining (symmetric) dissipation bracket in terms of the potential gradients. This last part is facilitated by using the constitutive equations provided by Eq. (37) to solve for the fluxes in terms of gradients of the potentials; however, the resulting relations are complicated for the general case. However, things are simpler in the case of the binary gas where the above relationships simplify considerably

leading to the following relations:

If we now use the suggested replacements by Eq. (55), Eqs. (51) (with the approximation $\alpha 12\u2248\alpha \u030212$) and (39) we get

This results in the following symmetric dissipation bracket:

The important item to notice from this exercise is that not only do we obtain a the required bilinear form for the dissipation bracket (e.g., Ref. 65) in terms of the gradient of the chemical potential differences and gradients of the temperature; but that all transport parameters involved are now explicitly expressed in terms of the true transport parameters of the problem (i.e., $D12,\u2009\lambda ,\u2009D1T,\u2009D2T$) with explicit expressions also provided by kinetic theory.^{12} Even if these correlations have been derived for binary mixtures before^{65} their correspondence is much clearer, and the approach illustrated here applies to more general, multidimensional systems. Furthermore, this new approach allows for added corrections that may be important in certain systems (such as that involved in the differentiation between the parameters $\alpha 12$ and $\alpha \u030212$). Moreover, this approach and the illustration here are valuable for extending the theory to binary systems applied beyond dilute gases and multicomponent cases, where it is anticipated that the unknown transport coefficients can much better fitted to experimental data using simple empirical models that modify the limiting results obtained from the kinetic theory of ideal gases, rather than starting from a general description that includes all dependences in terms of the state variables open for evaluation. We therefore believe that this connection to the kinetic theory of dilute gases is very valuable for future investigations of multicomponent heat and mass transfer in more complex systems.

## IV. BINARY DIFFUSION OF A DILUTE SOLUTE IN A LIQUID SOLVENT

As an illustrative flux-based transport application, it is instructive to consider a simple time-dependent, binary dilute liquid diffusion case, allowing for the time evolution of the full time-dependent solute species momenta. The governing equations in this case are: (1) the species mass densities evolution equation, Eq. (22), applied for species A (solute) assuming zero overall flow, which becomes

and (2) the corresponding species momentum equation, Eq. (60) for a stagnant medium, for which $v\xaf\u2248v\xafB=0\xaf\u21d2v\u2322\xafA=v\xafA$ and further, by neglecting higher order contributions, yields the following linearized equation:

If we now take the time derivative of the first equation, Eq. (65) and we subtract from it the gradient of the second equation, Eq. (66), we obtain a second order, linear wave equation

with characteristic speed

i.e., a correspondence is established in this limit between solute mass transfer and the wave equation with finite speed of propagation! For a typical low molecular weight liquid (e.g., water, $MA=18\u2009kg/kmol$) and for room temperature, $T=298\u2009K$

which is quite appreciable and of the same order of magnitude as the speed of sound in the medium (for water is equal to about 1500 m/s). Note that using this characteristic speed, which is comparable to the average Maxwell–Boltzmann root mean square (RMS) velocity $vrms=(3RGT)MA$, enables to rewrite Eq. (67) as a hyperbolic equation

which reduces to the usual, linear parabolic diffusion equation when $c2\u2192\u221e$. Therefore, with the new formalism we validate the traditional approach for most applications, except possibly those corresponding to very fast transients.

With this opportunity, we would like to clarify here that this is neither the first time to develop a hyperbolic diffusion equation to transport problems nor the above modeling procedure is the only way to accomplish it. As a matter of fact, there are many ways to obtain a hyperbolic diffusion equation in transport, the simplest one being adding a rate term to the traditional transport relations as first demonstrated for momentum in viscoelasticity by Maxwell^{66} and later in heat transfer by Cattaneo.^{67} In mass diffusion, this task has been so far mostly accomplished with the use of an internal variable—see, for example, Refs. 40 and 68 and references therein—although there have been interesting works to develop it from a microscopic point of view using Kramers equation for a Brownian particle.^{69} However, what we want to stress here is the universality of the approach through which this equation is obtained, applicable to all systems, with the only difference being in the values of the parameters involved. This restores physical expectations from the model, as for example that for a finite speed of signal propagation, without requiring us to impose them *a priori*. We feel therefore that this deserves a special discussion of some of the practical implications resulting from the use of such an equation, further illustrated below in a specific example of 1D transient mass transfer.

It is also worthwhile to explore Eq. (70), when applied to traditional diffusion problems to observe the time and length scales where the differences between species concentrations predicted by the parabolic and hyperbolic diffusion equations become prominent. It is also illustrative to review the solution methodologies that can be employed to solve the second order linear damped wave equation. Equation (70) can be rearranged to resemble the traditional parabolic diffusion equation or Fick's second law, albeit with a second order time derivative term as shown

Here, $\tau $ is the relaxation time scale for propagation defined as $\tau =DABc2$ that corresponds to the response time for diffusion when a concentration gradient is imposed. The above equation resembles the damped version of Fourier's law first developed by Cattaneo^{67} to account for inertia arising from finite speed of propagation in heat transfer.

This problem was further explored by Sieniutycz and Berry^{70} using variational principles of Onsager type. These variational principles subject the system to a constraint such that the growth of the total entropy generation is minimum. By working out vector couplings occurring in simultaneous heat and mass transfer, the authors were able to obtain a complete set of hyperbolic equations describing heat and mass transfer. In a previous work^{71} Sieniutycz presented relaxation times for several processes and systems. The relaxation times vary by several orders of magnitude depending on the length scale of the system; for instance, in Brownian diffusion, the relaxation time ranges from $3\xd710\u22128$ to $3\xd710\u22122$ s. As the relaxation time becomes smaller, so does the difference between predictions made by Fickian and hyperbolic diffusion equations.

For the purposes of illustration, consider a dilute solute species $A$ diffusing in a solvent $B$ in one-dimension through a slab of length $2L$, where the solute enters through both boundaries (applied as a step function) at $t=0$ with weight fraction $wA(1)$, with the solution initially being at solute weight fraction of $wA(0)$. Let us further consider as spatial domain of interest the interval $[\u22122L,0]$. For this system, the governing equation can be expressed by non-dimensionalizing Eq. (71) using the weight fraction at boundary $wA(1)$, the initial weight fraction $wA(0)$, the characteristic time $\tau $, and the corresponding diffusion length $DAB\tau $ as scales to construct appropriate dimensionless variables (denoted here with * superscripts), such that $x*=xDAB\tau $, $L*=LDAB\tau $, $wA*=wA\u2212wA(1)wA(0)\u2212wA(1)$ and $t*=t/\tau $.

The dimensionless system length $L*$ can then be considered as another measure to characterize the system. The differences between predictions of parabolic and hyperbolic equations are significant when the $L*\u2264O(1)$. Sandler and Dahler^{59} argue that the differences are only significant at length scales of $100\u2009\xc5$ in liquids or $105\u2009\xc5$ in gases.

Expressing the dimensional variables in Eq. (71) in terms of the above-defined dimensionless equivalents gives us

Exploring the symmetry over the centerline, this problem can be solved in half domain $(\u2212L*,0)$ and mirrored about $x*=\u2212L*$. We can therefore replace the left Dirichlet boundary condition at $x*=\u22122L*$ with a Neumann boundary condition at $x*=\u2212L*$ such that $\u2202wA*\u2202x*(x*=\u2212L*,t*)=0$. One additional boundary condition is required to obtain the coefficients, since the equation in second order in time, for which we propose $\u2202wA*\u2202t*(x*,t*=0)=0$. Solving the above equation using Fourier eigenfunctions such that $wA*=\u2211n=0\u221e\Theta n(t*)\Phi n(x*)$, the eigenfunctions and eigenvalues, $\Phi n(x*),\u2009\lambda n$, respectively, based on this combination of homogeneous Neumann and Dirichlet boundary conditions come out to be

where $n=0,1,2,\u2026$. Using a linear projection of variables along the eigenfunctions, defined as $\Theta n(t*)=\u222b0L*\Phi n(x*)wA*\u2009dx$, the transformed equation and boundary conditions become

The solution for this transformed equation is

where $\beta =1\u22124\lambda n2\tau *2$. Combining Eqs. (73) and (75), we can construct the solution for the dimensionless variable $wA*$ as follows:

which can be expressed in the dimensional space as

Figure 1 shows the above solution plotted for an initial solute weight fraction $wA(0)=0$ and weight fraction at the boundary $wA(1)=1$. At smaller times, one can observe the progression of the solute front across the spatial dimension. The dotted lines represent the solution of parabolic (Fickian) diffusion problem for the same boundary conditions for comparison. At short times, the parabolic diffusion equation has finite value everywhere, which is nonphysical since solute species cannot travel at infinite speed. This nonphysical behavior of the Fickian diffusion can be clearly contrasted with the improved inertial diffusion model that accounts for the finite speed of propagation in medium.

Another method to obtain an analytical solution for this problem valid at small times is by first considering the transient solution in an infinite domain ($\u2212\u221e\u2264x*\u2264\u221e$) of an initial delta function, i.e., subject to the initial conditions,

The solution to this problem has been developed by Gomez^{72} based on modified Bessel functions of the first kind. The solution was found to be

where $I0(x),\u2009I1(x)$ are the modified Bessel functions of the first kind and order 0, 1, respectively, defined using a series expansion as

Given the form of initial conditions used (78), the above solution can be considered as a Green's function, and the solution for any other initial value problem

can therefore be obtained as a convolution

For the semi-infinite domain ($x*\u22640$) and the Dirichlet boundary condition applied on the right boundary [$wA(x*=0,t*)=1$], the corresponding *f* function is

Furthermore, because for all values of $x*\u2264\u2212t*$ the solution is equal to 0, Eq. (82) can be simplified to

Evaluating the integral using Eq. (79) gives us the solution for the semi-infinite problem in the following form:

The integrals involved in the expression can be computed either numerically (using Gaussian quadrature) or analytically using truncated series expansion and integration of polynomial expressions—see Appendix C.

Yet another approach to solve the problem of hyperbolic diffusion equation is numerically by using finite differences. In this case, we first need to identify the characteristics of this problem along which the solution propagates, namely $x*=t*$ and $x*=\u2212t*$. The numerical solution of hyperbolic problems is restricted by the domain of dependence of the corresponding characteristics. In this case, stability requires (when explicit finite differences are used in time) the time steps ($\Delta t*$) to be smaller or equal than the spatial steps ($\Delta x*$). For simplicity and for illustration purposes, we propose here a first order finite difference scheme with the time steps ($\Delta t*$) equal to the spatial steps ($\Delta x*$). This can be most efficiently implemented using staggered grid points for discretization, as shown in Fig. 2. To properly apply the finite differences approximations, we first need to define the following material derivative operators:

where $D+$ is computed along the $x*=t*$(positive) characteristic, and $D\u2212$ is computed along the $x*=\u2212t*$(negative) characteristic. Equation (72) can then be written in terms of these operators as

Using the stencil shown in Fig. 2, the above expression can be discretized as follows:

Of course, when a higher accuracy is desired and/or the resolution of additional features present in the equations (additional terms, higher dimensions, and/or nonlinear dependencies) one needs to resort to more sophisticated numerical methods. Again, the equations being hyperbolic allows for a wide selection of really efficient numerical methods—see, for example, Refs. 37 and 73. A similar simplification is carried over to the specification and implementation of the boundary/initial conditions as discussed next.

Second, the propagation from the boundaries/initial conditions of any discontinuities needs to be evaluated numerically or analytically as appropriate for the problem at hand corresponding to the characteristic structure of the governing equation, Eq. (87). For example, in the case of the solution undergoing a step discontinuity from 0 to 1 at the point $x*=t*=0$, this step change propagates along the negative characteristic emanating from this point, satisfying the following equation:

which readily admits the analytical solution $\Delta w=exp\u2009(\u2212t2)$. After upgrading the values at the grid points due to the propagation of the boundary discontinuities, the discretized equation can then be iteratively solved for grid point $(i,j)$ recursively for each layer of time at a time.

Characteristic results are shown, in comparison to the available analytical solutions, in Figs. 3 and 4. In Fig. 3, a comparison of the solutions obtained from finite differences with the analytical solution for semi-infinite domain ($x*\u22640$) obtained in Eq. (85) for various spatial resolutions. A first order convergence of the numerical solution in terms of the common length/time step size is clearly demonstrated. The finite difference and analytical (series) solution to the hyperbolic diffusion problem in a slab has also been compared with numerical and analytical (series) solutions to the parabolic diffusion problem in Fig. 4. As seen there, the numerical solution to the hyperbolic problem shows comparable error to the numerical solution of the parabolic as compared to the common hyperbolic/parabolic solution achieved at times only one order of magnitude larger than the characteristic hyperbolic time $\tau $.

For transient problems involving very large length scales, the kind one would typically encounter in civil and environmental engineering, it becomes very critical to model the propagation front of diffusing species. Gomez *et al.*^{37} have shown how one can leverage the hyperbolic diffusion equation to take advantage of powerful numerical techniques like the discontinuous Galerkin (DG) approach to model convective and diffusive transport in systems involving such disparate length scales. The authors present an example of accidental spillage in a port, demonstrating the capability offered by hyperbolic equations for solving large scale problems efficiently with reduced computational effort.

The capability to express a parabolic problem in a hyperbolic form can therefore be significant even in cases where the hyperbolicity does not directly influence the results, i.e., when $L*\u226b1$, as it allows for the use of efficient hyperbolic numerical solvers. Of course, when $L*\u226b1$ and explicit solvers are used the stability requirements on the use of smaller time steps may be excessive. However, in that case, one can use an artificially high relaxation time to alleviate those restrictions, as long as its use does not affect appreciably the results. This choice is not different from the use of artificial diffusion or viscosity in complex fluid flow problems (like viscoelastic and/or turbulent^{74,75})

The linear hyperbolic diffusion equation, as explained before, is not new—it has appeared several times before in the literature been derived using a variety of approaches, ranging from kinetic theory,^{59,60} to variational approaches.^{70,71} What is new though is the systematic development of this equation within the context of nonequilibrium thermodynamics and the capability to extend it to accommodate nonlinear, flow, and/or multicomponent species effects.

## V. CONCLUSIONS

In this work, we demonstrated how the use of a systematic nonequilibrium thermodynamics approach, based on four simple steps (selection of state variables, description of the total energy (Hamiltonian), definition of the reversible (Poisson), and irreversible (dissipative) bracket dynamics), can systematically lead without *ad hoc* assumptions to the complete model equations of complex fluid flow, heat, and mass transfer, including the effects arising from the inertial dynamics of the mass fluxes, within a multicomponent system. It thus establishes to the fullest the analogy and unification between all transport processes, momentum, heat, and mass transfer, toward which R. B. Bird contributed so much with his work,^{21} making the present work an appropriate tribute to the celebration of his life. Simultaneously, we hopefully showed how his work could be amplified with judicious help from nonequilibrium thermodynamics. Although Bird was clearly aware and a good user of the linear irreversible thermodynamics [see, for example, the last chapter, Chap. 24, of Ref. 21 and an expert in the use of equilibrium thermodynamics to transport processes (as his last paper,^{15} easily demonstrates) curiously he had not exploited by himself much of the power of the reversible–irreversible structure of either SGBF-NET or GENERIC, albeit he did quote those works in the above-mentioned Chap. 24 of Ref. 21, most likely foreseeing their positive future impact to this field.

Several additional accomplishments are to be noted: First, we demonstrate how to correctly incorporate inertia (flux) effects in the modeling of coupled heat and mass transfer in a method consistent with nonequilibrium thermodynamics while preserving Galilean invariance. This approach incorporates species flux effects through the explicit use of their reduced momenta densities as additional variables. Thereby, this approach resolves a long-standing conundrum in nonequilibrium thermodynamics. Second, this modeling is shown to be successful in extending traditional descriptions of multicomponent diffusion–convection problems leading naturally to general expressions for diffusion and mass transfer without ambiguities or arbitrary assumptions. Those expressions correctly reduce to the extended (including heat transfer) Maxwell–Stefan equations and Fickian diffusion in the corresponding inertial-less limits of ideal gas and binary liquid with dilute solute, respectively. A successful comparison against kinetic theory predictions, in addition to providing for an independent validation of the modeling approach, enables obtaining explicit expressions for all the transport coefficients involved for this limiting behavior, while also providing a paradigm for extending the kinetic theory results to denser and multicomponent systems.

There are many open problems in multicomponent systems where the above-presented analysis may find important application. First, is the issue how to deal effectively with nonidealities and fugacities in extending the dilute gas kinetic theory relations for the transport coefficients. Moreover, as multicomponent systems present critical points where possibly molecular or thermal diffusion coefficients can vanish there is controversy how this can affect the transport properties and what are the limitations in observed behavior given the Onsager reciprocal relations.

Working directly with the mass fluxes (reduced momenta) and their dynamics has distinct advantages: On the one hand, this approach represents the physics more faithfully, allowing for finite speeds of propagation for the disturbances as well as for a wider variety of boundary conditions (described now in terms of the vector fluxes, not on scalar concentrations). On the other hand, it leads to fully hyperbolic equations that are easier to solve numerically through the use of powerful hyperbolic numerical solvers. The traditional parabolic convection–diffusion equation can be obtained in the limit of small relaxation times from its fully hyperbolic, finite signal propagation, equivalent. It is therefore proposed and illustrated here for a specific simple example, how to exploit the hyperbolic formulation resulting from the flux modeling of diffusion–convection problems in order to apply powerful numerical hyperbolic solvers to obtain solutions in their parabolic limits. Additional applications in many areas of chemical engineering are anticipated.

## ACKNOWLEDGMENTS

The authors acknowledge funding for this research provided by the National Science Foundation through CBET No. 1804911.

## AUTHOR DECLARATIONS

### Conflict of Interest

The authors have no conflicts to disclose.

## DATA AVAILABILITY

The data that support the findings of this study are available from the corresponding author upon reasonable request.

## NOMENCLATURE

### Description

- $aj$
$jth$ species/component in the system

- $c$
molar concentration

- $cs$
characteristic speed

- $C\xaf\xaf(k)$
conformation tensor corresponding to component $k$

- $dV$, $d3x$
differential volume

- $Dij$
molecular diffusion coefficient between $ith\u2212jth$ species

- $DiT$
Thermal diffusion coefficient of $ith$ species

- $h$
Hamiltonian functional density

- $H$
Hamiltonian functional

- $I\xaf\xaf$
identity tensor

- $j\xaf$
mass-based diffusive flux of $jth$ species

- $kT$
thermal diffusion ratio

- $L$
length

- $Mj$
molar mass/molecular weight

- $p$
thermodynamic pressure

- $q\xaf$
heat flux

- $Q\xaf\xaf$
viscous dissipation tensor

- $RG$
universal gas constant

- $s$
total entropy density

- $t$
time

- $T$
temperature of the system

- $u\xaf$
total momentum density

- $u\xaf(j)$
$jth$ species momentum density, $j\u22600$

- $u\u2322\xaf(j)$
$jth$ species reduced momentum density, $j\u22600$

- $u\u2322\xaf(0)$
heat momentum density

- $v\xaf$
velocity field of the system or reference frame

- $v\xaf(j)$
$jth$ species velocity, $j\u22600$

- $v\u2322\xaf(j)$
$jth$ species reduced velocity, $j\u22600$

- $v\u2322\xaf(0)$
heat velocity

- $wA$
weight fraction of species A

- $x$
spatial coordinate

- $xj$
molar fraction of $jth$ species

### Greek symbols

- $\alpha ij$
generalized coefficient of mass transfer between components $i$ and $j$, $i,j\u22600$

- $\alpha \u0302ij$
reduced mass transfer coefficient, $i,j\u22600$

- $\alpha i0$
Soret–Dufour coefficient for component $i$, $i\u22600$

- $\alpha 00$
generalized heat conductivity coefficient

- $\epsilon $
energy density functional

- $\kappa \u2032$
bulk viscosity

- $\lambda $
thermal conductivity

- $\mu $
shear viscosity

- $\mu j$
$jth$ species chemical potential

- $\mu 0j$
reference chemical potential of $jth$ species

- $\Xi $
primitive dissipation functional

- $\rho $
total system density

- $\rho j$
$jth$ species density

- $\rho 0$
effective heat mass density

- $\sigma \xaf\xaf$
stress tensor

- $\sigma s$
volumetric rate of entropy production

- $\tau $
characteristic time in hyperbolic diffusion

- $\Omega $
domain of interest in $x\xaf$ space

### Brackets/scripts

- $[\u22c5,\u22c5]$
dissipation bracket

- $[\u22c5,\u22c5]p$
primitive dissipation bracket

- $[\u22c5,\u22c5]ec$
entropy correction term

- $[\u22c5,\u22c5]a$
anti-symmetric part of dissipation bracket

- $[\u22c5,\u22c5]sp$
symmetric part of dissipation bracket

- ${\u22c5,\u22c5}$
Poisson bracket

- $\u27e8\u22c5\u27e9$
average/mean

- $\u22c5\u2009\Delta \xaf,\u2009\u22c5\u2009\xaf\xaf\Delta $
lower-convected Oldroyd time derivative for a tensor

- $\u2009\u22c5\u2009\u2192$
quantity in

*n*-vector space- $\u2009\u22c5\u2009\u2192\u2192$
matrix in $n\xd7n$ space

### APPENDIX A: CORRESPONDENCE BETWEEN THE OVERSCRIPT CURVED ARC AND REGULAR PRIMITIVE VARIABLES AND THEIR CORRESPONDING VOLTERRA DERIVATIVES OF AN ARBITRARY FUNCTIONAL *G*

For convenience and easy reference, we summarize here the correspondence between the two alternative $w\xaf$ and $w\u2322\xaf$ equivalent frameworks. First, it is clear that only the momentum density variables are affected. Second, from the definition of the overscript curved arc variables, Eqs. (6) and (5), their correspondence with the original momentum species densities is as follows:

Similarly, the correspondence between the original and overscript curved arc variables is

From those two sets of relationships, the following relationships can be derived between the corresponding Volterra derivatives of an arbitrary functional *G* using differentiation by parts:

and

### APPENDIX B: EQUIVALENCE BETWEEN THE MOST GENERAL AND THE PROPOSED BILINEAR DESCRIPTIONS FOR THE (PRIMITIVE) SYMMETRIC DISSIPATION BRACKET IN MULTICOMPONENT DIFFUSION

Considering the variables $w\u2322\xaf\u2261[\rho 1,\rho 2,\u2009\u2026\u2009,\rho n,u\xaf,u\u2322\xaf(1),u\u2322\xaf(2),\u2026\u2009,u\u2322\xaf(n),s]T$, the symmetries (parity) upon time reversal, and their physical interpretation [i.e., the definition of the reduced momenta based on the difference of the species velocity from the mass average velocity, Eq. (6)] the most general bilinear description for the (primitive) symmetric dissipation bracket for an isotropic medium is

where $Q\alpha \beta \gamma \epsilon \u2261\mu (\delta \alpha \gamma \delta \beta \epsilon +\delta \alpha \epsilon \delta \beta \gamma )+\kappa \u2032\delta \alpha \beta \delta \gamma \epsilon $ with $\mu ,\u2009\kappa \u2032$ are the shear and bulk viscosity^{41} and $aij,\u2009\kappa $ are the *i-j* generalized binary mass transfer coefficients and a scalar heat conductivity corresponding to this problem. With the caveat that when anisotropicity is present the $aij,\u2009\kappa $ need to be substituted by symmetric second order tensors, while $Q\alpha \beta \gamma \epsilon $ can have a more complex, anisotropic, structure, the above expression can easily be shown to present the most general bilinear description involving the above-mentioned variables. Of course, higher order dependencies on the Volterra derivatives of the Hamiltonian are also possible to appear in the dissipation bracket, still, close to equilibrium, upon linearization, those expressions also reduce to the above one.

At first sight, Eq. (B1) is different from the expression ultimately used, Eq. (15). We will show in this appendix though that they are equivalent when one considers that the reduced momenta are not independent but that they have to satisfy the constraint specified by Eq. (7). That constraint is easy to show that can also be described by the inner product of the *n*-vector $u\u2322\xaf\u2192\u2261[u\u2322\xaf(1),u\u2322\xaf(2),\u2009\u2026\u2009,u\u2322\xaf(n)]T$ with the unit *n*-vector $k\u2192\u22611n[1,1,\u2009\u2026\u2009,1\ufe38n\u2009terms]T$ where the Euclidean L2 norm is used for normalization

i.e., the *n*-vector $u\u2322\xaf\u2192\u2261[u\u2322\xaf(1),u\u2322\xaf(2),\u2009\u2026\u2009,u\u2322\xaf(n)]T$ is restricted to be in a subspace that is perpendicular to with the unit *n*-vector $k\u2192\u22611n[1,1,\u2009\u2026\u2009,1\ufe38n\u2009terms]T$. Due to this constraint, it is well known that a similar condition also applies^{41} to the *n*-vector of the Volterra derivatives with respect to the reduced momenta of any arbitrary functional *G* (G = *H*, *F*, ….), $\delta G\delta u\u2322\xaf\u2192\u2261[\delta G\delta u\u2322\xaf(1),\delta G\delta u\u2322\xaf(2),\u2009\u2026\u2009,\delta G\delta u\u2322\xaf(n)]T$

i.e., the *n*-vector $\delta G\delta u\u2322\xaf\u2192\u2261[\delta G\delta u\u2322\xaf(1),\delta G\delta u\u2322\xaf(2),\u2009\u2026\u2009,\delta G\delta u\u2322\xaf(n)]T$ is restricted to be in a subspace that is perpendicular to the unit *n*-vector $k\u2192\u22611n[1,1,\u2009\u2026\u2009,1\ufe38n\u2009terms]T$.

Due to this constraint, not all of the *n*(*n *+* *1)/2 components of (a symmetric, due to Onsager/Casimir reciprocity relations^{4–6}) *n* × *n* coefficient matrix $a\u2192\u2192$ are independent. Instead, the only relevant information (that can lead to measurable results) belongs to the subspace of the *n* × *n* coefficient symmetric matrix $a\u2192\u2192$ that satisfy the corresponding to Eq. (B3) condition

i.e., without loss of generality we can assume that the *n* × *n* coefficient symmetric matrix $a\u2192\u2192$ has a zero eigenvalue along the unit *n*-vector $k\u2192\u22611n[1,1,\u2009\u2026\u2009,1\ufe38n\u2009terms]T$ which becomes one of its eigenvectors.

The most general expression that satisfy this condition can be shown to be provided (through a partial eigenvalue analysis, always possible for a symmetric matrix) by

where $[q\u2192(1),\u2009q\u2192(2),\u2026,q\u2192(n\u22121)]$ is a set of *n*-1 independent orthonormal *n*-vectors perpendicular to the unit *n*-vector $k\u2192\u22611n[1,1,\u2009\u2026\u2009,1\ufe38n\u2009terms]T$, $0\u2322\u2192$ the n-1 zero vector and $a\u2322\u2192\u2192$ a general, symmetric $(n\u22121)\xd7(n\u22121)$ matrix. It can easily be confirmed that $a\u2192\u2192$ constructed as indicated by Eq. (B5) duly satisfies the conditions of Eq. (B4).

Although Eq. (B5) is helpful in identifying one way to construct $a\u2192\u2192$ as well as in constructively showing that the actual space where $a\u2192\u2192$ belongs has *n*(*n*-1)/2 dimensions [the number of independent components of the general, symmetric $(n\u22121)\xd7(n\u22121)$ matrix $a\u2322\u2192\u2192$] it does not offer neither the easiest nor the most physical way to define it as its components by necessity depend on the selection of the orthonormal vectors $[q\u2192(1),\u2009q\u2192(2),\u2026,q\u2192(n\u22121)]$. However, once we have identified, thanks to Eq. (B5) the dimensionality of the space of $a\u2192\u2192$ suffices in order to uniquely define it to decompose it in terms of any set of *n*(*n*-1)/2 independent members of that space. As it turns out, those are most naturally constructed based on the decomposition offered in Eq. (15) using the corresponding to each $\alpha ij$ coefficients $(i=1,2,\u2026,n,\u2009j=1,2,\u2026,n,\u2009j\u2260i)$ *n* × *n* elementary matrices $a\u2192\u2192ij$ defined through the specification of their k, l component, $(a\u2192\u2192ij)kl$, $(k=1,2,\u2026,n,\u2009l=1,2,\u2026,n)$ as

It is indeed straightforward to show that each one of those matrices satisfies the conditions corresponding to Eq. (B4) (i.e., it belongs to the right space) and that for $(i,j)\u2260(j,i)$ they are independent, and therefore, as they are *n*(*n*-1)/2 in number they form a basis for that space in terms of which the generic matrix $a\u2192\u2192$ can be described, as

### APPENDIX C: THE INTEGRALS APPEARING EQ. (85) CAN BE SOLVED ANALYTICALLY BY FIRST PERFORMING SERIES EXPANSION OF THE MODIFIED BESSEL FUNCTIONS OF THE FIRST KIND, FOLLOWED BY INTEGRATION OF THE POLYNOMIAL TERMS THAT APPEAR IN THE SERIES EXPANSION

For the first integral term $\u222b\u2212txI0(z\u0302)\u2009dx\u0302$, the expansion can be expressed as follows:

Similarly, the second integral term can be expressed as

The above integrals can be used to express Eq. (85) as a series expansion in the following manner:

## References

*:*