Flux-based modeling of heat and mass transfer in multicomponent systems

Beris, Antony N.; Jariwala, Soham; Wagner, Norman J.

doi:10.1063/5.0085444

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).² 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)²,² 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¹¹ and was finalized shortly afterward in the monumental treatise of “molecular theory of gases and liquids” by Hirschfelder, Curtiss, and Bird.¹² 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¹³—in a recent article by Bird and Klingenberg¹⁴ 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,¹⁵ 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¹⁶ 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.¹⁷ 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.¹⁷ 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.¹⁷ 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¹⁸ 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.³² 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³⁶) 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,³⁷ 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).⁴¹ This represents one of the simplest formulations of the most modern implementations of nonequilibrium thermodynamics, the GENERIC formulation⁴² to which it has been shown to be completely equivalent when applied to macroscopic modeling.⁴³ 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.⁴⁴ 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.¹² 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.¹² 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¹² 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.⁴² 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⁴⁶ that outlines its influence in elastic fluid modeling—for a more detailed presentation the interested reader is referred to the above-mentioned research monograph.⁴¹

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⁴⁹ and Arnold.⁵⁰ 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, ${\begin{matrix} z_{1} & z_{2} & \dots & z_{n} \end{matrix}}$ ⁠, 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 = \int_{V} h (\begin{matrix} z_{1} & z_{2} & \dots & z_{n} \end{matrix}; \nabla z_{i}; \nabla^{2} z_{i}; etc .) d V$ ⁠, 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.⁴¹
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.⁴¹ 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

\frac{d F}{d t} = {F, H} + [F, H] .

(1)

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

\frac{d F}{d t} = \int_{V} \sum_{i = 1}^{n} \frac{δ F}{δ z_{i}} \frac{\partial z_{i}}{\partial t} d V,

(2)

and equating the factors in front of each Volterra derivative $\frac{δ F}{δ z_{i}}$ ⁠.⁴¹ 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, $a_{j}$ ⁠, where $j = 1, \dots, n$ ⁠, the standard set of variables are: (a) each one of the n species mass densities, $ρ_{j}$ ⁠, and (b) their corresponding momentum densities

{\underline{u}}^{(j)} \equiv ρ_{j} {\underline{v}}^{(j)}, j = 1, \dots, n,

(3)

where ${\underline{v}}^{(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

\underline{u} \equiv ρ \underline{v} \equiv \sum_{j = 1}^{n} ρ_{j} {\underline{v}}^{(j)},

(4)

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

ρ = \sum_{j = 1}^{n} ρ_{j}; \underline{u} = \sum_{j = 1}^{n} {\underline{u}}^{(j)} .

(5)

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 $ρ_{j}$ [through which the total mass density ρ, can also be obtained, from the first relationship in Eq. (5)] (b) the total momentum density, $\underline{u}$ ⁠, (c), the reduced momentum densities corresponding to each species

{\underline{\overset{⌢}{u}}}^{(j)} \equiv ρ_{j} ({\underline{v}}^{(j)} - \underline{v}) \equiv ρ_{j} {\underline{\overset{⌢}{v}}}^{(j)},

(6)

(d) the total entropy density, $s$ ⁠. Note that in this case, the total momentum density (⁠ $\underline{u}$ ⁠) 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:

\sum_{j = 1}^{n} {\underline{\overset{⌢}{u}}}^{(j)} = \underline{0} .

(7)

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 ${\underline{\overset{⌢}{u}}}^{(0)}$ ⁠. This is defined in terms of the heat flux $\underline{q}$ as

{\underline{\overset{⌢}{u}}}^{(0)} \equiv ρ_{0} {\underline{\overset{⌢}{v}}}^{(0)} \equiv \frac{ρ_{0}}{s T} \underline{q},

(8)

where the effective heat mass density, $ρ_{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 ${\underline{\overset{⌢}{v}}}^{(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 ${\underline{\overset{⌢}{v}}}^{(i)}$ does with respect to the diffusive flux of species i. After the addition of this new reduced heat momentum density, the appropriate state $2 n + 3$ variable space can be defined as $\underline{\overset{⌢}{w}} \equiv {[s, ρ_{1}, \dots, ρ_{n}, \underline{u}, {\underline{\overset{⌢}{u}}}^{(0)}, {\underline{\overset{⌢}{u}}}^{(1)}, \dots, {\underline{\overset{⌢}{u}}}^{(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 $2 n + 2$ ⁠.

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

\begin{matrix} H = \int_{Ω} (\sum_{j = 1}^{n} \frac{1}{2} ρ_{j} {(v^{(j)})}^{2} + \frac{1}{2} ρ_{0} {({\overset{⌢}{v}}^{(0)})}^{2} + ε) d V \\ = \int_{Ω} (\frac{1}{2} ρ v^{2} + \sum_{j = 1}^{n} \frac{1}{2} ρ_{j} {({\underline{v}}^{(j)} - \underline{v})}^{2} + \frac{1}{2} ρ_{0} {({\overset{⌢}{v}}^{(0)})}^{2} + ε) d V \\ = \int_{Ω} (\frac{1}{2} \frac{u^{2}}{ρ} + \sum_{j = 1}^{n} \frac{1}{2} \frac{{({\overset{⌢}{u}}^{(j)})}^{2}}{ρ_{j}} + \frac{1}{2} \frac{{({\overset{⌢}{u}}^{(0)})}^{2}}{ρ_{0}} + ε) d^{3} x \\ = \int_{Ω} (\frac{1}{2} \frac{u^{2}}{\sum_{j = 1}^{n} ρ_{j}} + \sum_{j = 0}^{n} \frac{1}{2} \frac{{({\overset{⌢}{u}}^{(j)})}^{2}}{ρ_{j}} + ε) d^{3} x, \end{matrix}

(9)

where the last term is written in terms of the prescribed system variables with $ε = ε (s, ρ_{1}, \dots, ρ_{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, $ρ_{0} = ρ_{0} (s, ρ_{1}, \dots, ρ_{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, $c_{s}$ (in an equation, that is, introduced for convenience and shall be used as a proper definition for the characteristic speed, $c_{s}$ ⁠) as

ρ_{0} \equiv \frac{s T}{{(c_{s})}^{2}} .

(10)

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.⁵⁴ with the coefficient $z = \frac{q_{0}}{{(s T)}^{2}}$ representing the magnitude of an isotropic tensor $\underline{\underline{Z}} = z \underline{\underline{I}}$ 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

\begin{matrix} \frac{δ H}{δ ρ_{j}} = - \frac{1}{2} \frac{u^{2}}{ρ^{2}} - \frac{1}{2} \frac{{({\overset{⌢}{u}}^{(j)})}^{2}}{ρ_{j}^{2}} - \frac{1}{2} \frac{{({\overset{⌢}{u}}^{(0)})}^{2}}{ρ_{0}^{2}} \frac{\partial ρ_{0}}{\partial ρ_{j}} + \frac{\partial ε}{\partial ρ_{j}} = - \frac{1}{2} v^{2} - \frac{1}{2} {({\underline{\overset{⌢}{v}}}^{(j)})}^{2} - \frac{1}{2} {({\underline{\overset{⌢}{v}}}^{(0)})}^{2} \frac{\partial ρ_{0}}{\partial ρ_{j}} + \frac{1}{M_{j}} μ_{j}, \\ \frac{δ H}{δ \underline{u}} = \underline{v}, \\ \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(j)}} = ({\underline{v}}^{(j)} - \underline{v}) - \frac{1}{n} (\sum_{k = 1}^{n} {\underline{v}}^{(k)}) + \underline{v} = {\underline{v}}^{(j)} - \frac{1}{n} (\sum_{k = 1}^{n} {\underline{v}}^{(k)}), \\ \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(0)}} = \frac{1}{ρ_{0}} {\underline{\overset{⌢}{u}}}^{(0)} \equiv {\underline{\overset{⌢}{v}}}^{(0)} \equiv \frac{1}{s T} \underline{q}, \\ \frac{δ H}{δ s} \equiv T = - \frac{1}{2} {({\underline{\overset{⌢}{v}}}^{(0)})}^{2} \frac{\partial ρ_{0}}{\partial s} + \frac{\partial ε}{\partial s}, \end{matrix}

(11)

where, as usual, the jth chemical potential and temperature are defined as $μ_{j} \equiv M_{j} \frac{\partial ε}{\partial ρ_{j}}$ and j takes the values $j = 1, 2, \dots, 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 ${[\begin{matrix} 1 & 1 & \dots & 1 \end{matrix}]}^{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_{α β}^{(k)} = ⟨ ρ_{k} R_{α}^{(k)} R_{β}^{(k)} ⟩ \equiv ρ_{k} c_{α β}^{(k)}$ as the mass-weighted average of the second moment $c_{α β}^{(k)}$ of a structural position vector, $R_{α}^{(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, $h_{e} = h_{e} (ρ_{k}, {\underline{\underline{C}}}^{(k)}, T)$ that also contributes to the Volterra derivative $\frac{δ H}{δ ρ_{k}}$ an extra elastic energy related term $\frac{\partial h_{e}}{\partial ρ_{k}}$ ⁠. This effect can lead to very interesting phenomena, such as stress induced migration,²⁸ 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, $\underline{v}$ ⁠, corresponding to the total momentum, $\underline{u}$ ⁠. This is in contrast to the multi-fluid model.³⁵ The reduced momenta, ${\underline{\overset{⌢}{u}}}^{(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

\begin{matrix} {F, H} = - \sum_{i = 1}^{n} \int_{Ω} [\frac{δ F}{δ ρ_{i}} \nabla_{β} (\frac{δ H}{δ u_{β}} ρ_{i}) - \frac{δ H}{δ ρ_{i}} \nabla_{β} (\frac{δ F}{δ u_{β}} ρ_{i})] d^{3} x \\ - \sum_{i = 0}^{n} \int_{Ω} [\frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(i)}} \nabla_{β} (\frac{δ H}{δ u_{β}} {\overset{⌢}{u}}_{α}^{(i)}) - \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}} \nabla_{β} (\frac{δ F}{δ u_{β}} {\overset{⌢}{u}}_{α}^{(i)})] d^{3} x \\ - \sum_{i = 0}^{n} \int_{Ω} [{\overset{⌢}{u}}_{β}^{(i)} (\frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(i)}} \nabla_{β} (\frac{δ H}{δ u_{α}}) - \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}} \nabla_{β} (\frac{δ F}{δ u_{α}}))] d^{3} x \\ - \int_{Ω} [\frac{δ F}{δ u_{α}} \nabla_{β} (\frac{δ H}{δ u_{β}} u_{α}) - \frac{δ H}{δ u_{α}} \nabla_{β} (\frac{δ F}{δ u_{β}} u_{α})] d^{3} x \\ - \int_{Ω} [\frac{δ F}{δ s} \nabla_{β} (\frac{δ H}{δ u_{β}} s) - \frac{δ H}{δ s} \nabla_{β} (\frac{δ F}{δ u_{β}} s)] d^{3} x, \end{matrix},

(12)

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]}_{s p}$

{[F, H]}_{p} = {[F, H]}_{a} + {[F, H]}_{s p} .

(13)

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 $\underline{u}$ being replaced by the species velocities ${\underline{u}}^{(i)}$ and then the Volterra derivatives with respect to the species velocities being expressed in terms of the appropriate elements of the variable space $\underline{\overset{⌢}{w}}$ ⁠, i.e., the reduced species momenta ${\underline{\overset{⌢}{u}}}^{(i)}$ and the mass average velocity $\underline{u}$ [see Appendix A, Eq. (A3)].

The contributions of the Volterra derivatives with respect to the mass average velocity $\underline{u}$ 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

\begin{matrix} {[F, H]}_{a} = - \int_{Ω} \sum_{i = 1}^{n} (\begin{array}{l} \frac{δ F}{δ ρ_{i}} \nabla_{α} ((ρ_{i} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}}) - \frac{ρ_{i}}{ρ} (\sum_{k = 1}^{n} ρ_{k} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(k)}})) \\ - \frac{δ H}{δ ρ_{i}} \nabla_{α} ((ρ_{i} \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(i)}}) - \frac{ρ_{i}}{ρ} (\sum_{k = 1}^{n} ρ_{k} \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(k)}})) \end{array}) d^{3} x \\ - \int_{Ω} \sum_{i = 1}^{n} (\begin{array}{l} \frac{δ F}{δ s} \nabla_{α} ((s_{i} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}}) - \frac{s_{i}}{s} (\sum_{k = 1}^{n} s_{k} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(k)}})) \\ - \frac{δ H}{δ s} \nabla_{α} ((s_{i} \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(i)}}) - \frac{s_{i}}{s} (\sum_{k = 1}^{n} s_{k} \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(k)}})) \end{array}) d^{3} x - \int_{Ω} (\begin{array}{l} \frac{δ F}{δ s} \nabla_{α} (s \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}}) \\ - \frac{δ H}{δ s} \nabla_{α} (s \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(0)}}) \end{array}) d^{3} x . \end{matrix} .

(14)

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

\begin{matrix} {[F, H]}_{s p} = - \int_{Ω} \sum_{i, j = 1}^{n} (α_{i j} (\frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(i)}} - \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(j)}}) (\frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}} - \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}})) d^{3} x \\ - \int_{Ω} \sum_{i}^{n} α_{i 0} (\frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(i)}} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}} + \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}} \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(0)}}) d^{3} x - \int_{Ω} α_{00} \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(0)}} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}} d^{3} x \\ - \int_{Ω} Q_{α β γ ε} \nabla_{α} \frac{δ F}{δ u_{β}} \nabla_{γ} \frac{δ H}{δ u_{ε}} d^{3} x, \end{matrix}

(15)

where $Q_{α β γ ε} \equiv μ (δ_{α γ} δ_{β ε} + δ_{α ε} δ_{β γ}) + κ^{'} δ_{α β} δ_{γ ε}$ with $μ, κ^{'}$ are the shear and bulk viscosity⁴¹ and $α_{i j}, α_{i 0}, α_{00}$ are the i-j generalized binary mass transfer coefficients, $i, j = 1, 2, \dots n, j \neq i$ ⁠, ith generalized Soret/Dufour coefficients, $i = 1, 2, \dots, 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 $α_{i j}$ 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

\sum_{i = 1}^{n} α_{i 0} = 0.

(16)

In addition, it is important to note that as a result of the symmetry of the transport coefficients $α_{i j}$ 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 $α_{i j}, α_{i 0}, α_{00}$ ⁠; suffice to mention though here that in the particular case of a binary mixture those requirements reduce to (a) the transport coefficients $α_{12}, α_{00}$ to be non-negative and (b) the Soret/Dufour coefficients $α_{10} = - α_{20}$ to satisfy the inequality ${(α_{10})}^{2} \leq 2 α_{12} α_{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.⁵⁶

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.⁴¹ In general, if the primitive dissipative bracket is offered as

{[F, H]}_{p} \equiv \int_{Ω} Ξ (F, H) d^{3} x,

(17)

then, the entropy correction term, ${[F, H]}_{e c}$ is defined as

{[F, H]}_{e c} \equiv - \int_{Ω} \frac{1}{\frac{δ H}{δ s}} \frac{δ F}{δ s} Ξ (H, H) d^{3} x,

(18)

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

[F, H] \equiv {[F, H]}_{p} + {[F, H]}_{e c} .

(19)

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

\frac{d F}{d t} = \int_{Ω} \sum_{j} \frac{δ F}{δ {\overset{⌢}{w}}_{j}} \frac{\partial {\overset{⌢}{w}}_{j}}{\partial t} d^{3} x = \int_{Ω} (\sum_{i = 1}^{n} \frac{δ F}{δ ρ_{i}} \frac{\partial ρ_{i}}{\partial t} + \frac{δ F}{δ u_{α}} \frac{\partial u_{α}}{\partial t} + \sum_{i = 0}^{n} \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(i)}} \frac{\partial {\overset{⌢}{u}}_{α}^{(i)}}{\partial t} + \frac{δ F}{δ s} \frac{\partial s}{\partial t}) d^{3} x .

(20)

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, $\frac{d ρ_{i}}{d t}$ ⁠, would be equal to the coefficient of $\frac{δ F}{δ ρ_{i}}$ in the one generator expression for $\frac{d F}{d t}$ ⁠, based on the evaluation of the Volterra derivatives of the Hamiltonian provided by Eq. (11), leading to

\begin{matrix} \frac{\partial ρ_{i}}{\partial t} = \underset{From Poisson bracket}{\underset{︸}{- \nabla_{β} (\frac{δ H}{δ u_{β}} ρ_{i})}} \underset{From dissipation bracket}{\underset{︸}{- \nabla_{β} (\frac{δ H}{δ {\overset{⌢}{u}}_{β}^{(i)}} ρ_{i}) + \nabla_{β} \frac{ρ_{i}}{ρ} (\sum_{k = 1}^{n} \frac{δ H}{δ {\overset{⌢}{u}}_{β}^{(k)}} ρ_{k})}} \\ = - \nabla_{β} (v_{β} ρ_{i}) - \nabla_{β} ((v_{β}^{(i)} - \frac{1}{n} (\sum_{j = 1}^{n} v_{β}^{(j)})) ρ_{i}) + \nabla_{β} (\frac{ρ_{i}}{ρ} \sum_{k = 1}^{n} ((v_{β}^{(k)} - \frac{1}{n} (\sum_{j = 1}^{n} v_{β}^{(j)})) ρ_{k})) \\ = - \nabla_{β} (v_{β} ρ_{i}) - \nabla_{β} ((v_{β}^{(i)} - \frac{1}{n} (\sum_{j = 1}^{n} v_{β}^{(j)})) ρ_{i}) + \nabla_{β} (\frac{ρ_{i}}{ρ} \sum_{k = 1}^{n} v_{β}^{(k)} ρ_{k}) - \nabla_{β} (\frac{ρ_{i}}{ρ} (- \frac{1}{n} (\sum_{j = 1}^{n} v_{β}^{(j)})) \sum_{k = 1}^{n} ρ_{k}) \\ = - \nabla_{β} (v_{β} ρ_{i}) - \nabla_{β} ((v_{β}^{(i)} - \frac{1}{n} (\sum_{j = 1}^{n} v_{β}^{(j)})) ρ_{i}) + \nabla_{β} (\frac{ρ_{i}}{ρ} (ρ v_{β})) - \nabla_{β} (\frac{ρ_{i}}{ρ} (- \frac{1}{n} (\sum_{j = 1}^{n} v_{β}^{(j)})) ρ) \\ = - \nabla_{β} (v_{β} ρ_{i}) - \nabla_{β} ((v_{β}^{(i)} - v_{β}) ρ_{i}) = - \nabla_{β} (v_{β}^{(i)} ρ_{i}), \end{matrix}

(21)

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

\frac{\partial ρ_{i}}{\partial t} = - \nabla_{β} (v_{β}^{(i)} ρ_{i}) .

(22)

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

\begin{matrix} \sum_{i = 1}^{n} \frac{\partial ρ_{i}}{\partial t} = - \sum_{i = 1}^{n} \nabla_{β} (v_{β}^{(i)} ρ_{i}) \Leftrightarrow \frac{\partial ρ}{\partial t} = - \nabla_{β} (v_{β} ρ), \end{matrix}

(23)

where $ρ = \sum_{i = 1}^{n} ρ_{i}$ ⁠.

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

\frac{\partial u_{α}}{\partial t} = - \nabla_{β} (u_{α} v_{β}) + \nabla_{β} σ_{β α}^{(c)} - \nabla_{α} p + \nabla_{β} (Q_{α β γ ε} \nabla_{γ} v_{ε}),

(24)

where $σ_{β α}^{(c)}$ is a reversible contribution to the stress, as it arises from the reduced species momenta

σ_{β α}^{(c)} = - \sum_{j = 0}^{n} {\overset{⌢}{u}}_{β}^{(j)} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}} = - ρ_{0} \frac{q_{β} q_{α}}{{(s T)}^{2}} - \sum_{j = 1}^{n} ρ_{j} {\overset{⌢}{v}}_{β}^{(j)} {\overset{⌢}{v}}_{α}^{(j)},

(25)

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

p \equiv - h + s \frac{δ H}{δ s} + \sum_{i = 1}^{n} ρ_{i} \frac{δ H}{δ ρ_{i}} + u_{β} \frac{δ H}{δ u_{β}} + \sum_{i = 0}^{n} {\overset{⌢}{u}}_{β}^{(i)} \frac{\partial H}{\partial {\overset{⌢}{u}}_{β}^{(i)}} .

(26)

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

{p |}_{v - contributions} \equiv - \frac{1}{2} ρ v^{2} + \sum_{i = 1}^{n} ρ_{i} (- \frac{1}{2} \frac{u^{2}}{ρ^{2}}) + u_{β} (\frac{u_{β}}{ρ}) = ρ v^{2} (- \frac{1}{2} - \frac{1}{2} + 1) = 0.

(27)

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:

\frac{\partial v_{α}}{\partial t} = - v_{β} \nabla_{β} (v_{α}) + \frac{1}{ρ} (\nabla_{β} σ_{β α}^{(c)} - \nabla_{α} p + \nabla_{β} (Q_{α β γ ε} \nabla_{γ} v_{ε})) .

(28)

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

\frac{\partial s}{\partial t} = - \nabla_{β} (v_{β} s) - \sum_{i = 1}^{n} \nabla_{β} ({\overset{⌢}{v}}_{β}^{(i)} s_{i}) - \nabla_{β} ({\overset{⌢}{v}}_{β}^{(0)} s) + σ_{s},

(29)

where $σ_{s}$ represents the total rate of entropy production and is given as

σ_{s} = \frac{1}{T} (\sum_{i, j = 1}^{n} (α_{i j} {(v_{α}^{(i)} - v_{α}^{(j)})}^{2}) + 2 v_{α}^{(0)} (\sum_{i}^{n} α_{i 0} (v_{α}^{(i)} - \frac{1}{n} (\sum_{k = 1}^{n} v_{α}^{(k)}))) + α_{00} {(v_{α}^{(0)})}^{2} + \underline{\underline{\underline{\underline{Q}}}} : : \underline{\underline{\nabla v}} \underline{\underline{\nabla v}}) .

(30)

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:

\begin{matrix} \frac{\partial {\overset{⌢}{u}}_{α}^{(0)}}{\partial t} = - \nabla_{β} ({\overset{⌢}{u}}_{α}^{(0)} v_{β}) - {\overset{⌢}{u}}_{β}^{(0)} \nabla_{β} v_{α} - s \nabla_{α} \frac{δ H}{δ s} - \sum_{j = 1}^{n} (α_{j 0} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}}) - α_{00} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}}, \\ \frac{\partial {\overset{⌢}{u}}_{α}^{(i)}}{\partial t} = {\begin{array}{l} \begin{array}{l} - \nabla_{β} ({\overset{⌢}{u}}_{α}^{(i)} v_{β}) - {\overset{⌢}{u}}_{β}^{(i)} \nabla_{β} v_{α} - ρ_{i} \nabla_{α} \frac{δ H}{δ ρ_{i}} \\ + \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} \frac{δ H}{δ ρ_{j}}, \\ - α_{i 0} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}} - \sum_{j = 1}^{n} (α_{i j} + α_{j i}) (\frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}} - \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}}), \end{array} & i = 1, 2, \dots, n . \end{array} \end{matrix}

(31)

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

\begin{matrix} \overset{Δ}{{\overset{⌢}{u}}_{α}^{(0)}} = - s \nabla_{α} \frac{δ H}{δ s} - \sum_{j = 1}^{n} (α_{j 0} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}}) - α_{00} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}}, \\ \overset{Δ}{{\overset{⌢}{u}}_{α}^{(i)}} = - \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} (\frac{δ H}{δ ρ_{i}} - \frac{δ H}{δ ρ_{j}}) - α_{i 0} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}} - \sum_{j = 1}^{n} (α_{i j} + α_{j i}) (\frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}} - \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}}), i = 1, 2, \dots, n \end{matrix}

(32)

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

\overset{Δ}{w_{α}} \equiv \frac{\partial w_{α}}{\partial t} + \nabla_{β} (v_{β} w_{α}) + w_{β} \nabla_{β} (v_{α}) .

(33)

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

\begin{matrix} \overset{Δ}{{\overset{⌢}{u}}_{α}^{(0)}} = - s \nabla_{α} T - \sum_{j = 1}^{n} (α_{j 0} (v_{α}^{(j)} - \frac{1}{n} (\sum_{k = 1}^{n} v_{α}^{(k)}))) - α_{00} \frac{q_{α}}{s T}, \\ \overset{Δ}{{\overset{⌢}{u}}_{α}^{(i)}} = \begin{matrix} \begin{matrix} \begin{array}{l} - \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} (- \frac{1}{2} ({({\overset{⌢}{v}}^{(i)})}^{2} - {({\overset{⌢}{v}}^{(j)})}^{2}) - \frac{1}{2} {({\overset{⌢}{v}}^{(0)})}^{2} (\frac{\partial ρ_{0}}{\partial ρ_{i}} - \frac{\partial ρ_{0}}{\partial ρ_{j}}) + \frac{μ_{i}}{M_{i}} - \frac{μ_{j}}{M_{j}}) \\ - α_{i 0} \frac{q_{α}}{s T} - \sum_{j = 1}^{n} (α_{i j} + α_{j i}) (v_{α}^{(i)} - v_{α}^{(j)}), \end{array} \end{matrix} & i = 1, 2, \dots, n \end{matrix} . \end{matrix}

(34)

Note that at equilibrium (i.e., in the absence of any flow and inhomogeneities) the values of the reduced momenta densities $\underline{{\overset{⌢}{u}}_{i}}$ ⁠, and their lower convected time derivatives, are equal to $\underline{0}$ ⁠.

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)

\begin{matrix} \frac{\partial u_{α}^{(i)}}{\partial t} = - \nabla_{β} ({\overset{⌢}{u}}_{α}^{(i)} v_{β}) - {\overset{⌢}{u}}_{β}^{(i)} \nabla_{β} v_{α} - v_{α} \nabla_{β} (u_{β}^{(i)}) + ρ_{i} (- v_{β} \nabla_{β} (v_{α}) + \frac{1}{ρ} (\nabla_{β} σ_{β α}^{(c)} - \nabla_{α} p + \nabla_{β} (Q_{α β γ ε} \nabla_{γ} v_{ε}))) \\ - ρ_{i} \nabla_{α} \frac{δ H}{δ ρ_{i}} + \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} \frac{δ H}{δ ρ_{j}} - α_{i 0} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}} - \sum_{j = 1}^{n} (α_{i j} + α_{j i}) (\frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}} - \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}}) \\ = - \nabla_{β} (ρ {}_{i}{({\overset{⌢}{v}}_{α}^{(i)} v_{β} + v_{α} v_{β}^{(i)})}) + \frac{ρ_{i}}{ρ} (\nabla_{β} σ_{β α}^{(c)} - \nabla_{α} p + \nabla_{β} (Q_{α β γ ε} \nabla_{γ} v_{ε})) \\ - ρ_{i} \nabla_{α} \frac{δ H}{δ ρ_{i}} + \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} \frac{δ H}{δ ρ_{j}} - α_{i 0} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}} - \sum_{j = 1}^{n} (α_{i j} + α_{j i}) (\frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}} - \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}}) \\ = - \nabla_{β} (ρ {}_{i}{(v_{α}^{(i)} v_{β}^{(i)} - {\overset{⌢}{v}}_{α}^{(i)} {\overset{⌢}{v}}_{β}^{(i)})}) + \frac{ρ_{i}}{ρ} (\nabla_{β} σ_{β α}^{(c)} - \nabla_{α} p + \nabla_{β} (Q_{α β γ ε} \nabla_{γ} v_{ε})) \\ - ρ_{i} \nabla_{α} \frac{δ H}{δ ρ_{i}} + \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} \frac{δ H}{δ ρ_{j}} - α_{i 0} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(0)}} - \sum_{j = 1}^{n} (α_{i j} + α_{j i}) (\frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(i)}} - \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}}) . \end{matrix} .

(35)

The last equality was obtained by using the identities

ρ {}_{i}{({\overset{⌢}{v}}_{α}^{(i)} v_{β} + v_{α} v_{β}^{(i)})} = {\overset{⌢}{u}}_{α}^{(i)} v_{β} + (u_{α}^{(i)} - {\overset{⌢}{u}}_{α}^{(i)}) v_{β}^{(i)} = u_{α}^{(i)} v_{β}^{(i)} - {\overset{⌢}{u}}_{α}^{(i)} (v_{β}^{(i)} - v_{β}) = u_{α}^{(i)} v_{β}^{(i)} - {\overset{⌢}{u}}_{α}^{(i)} {\overset{⌢}{v}}_{β}^{(i)} = ρ {}_{i}{(v_{α}^{(i)} v_{β}^{(i)} - {\overset{⌢}{v}}_{α}^{(i)} {\overset{⌢}{v}}_{β}^{(i)})} .

(36)

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,¹⁷ the only difference is that we include a term $\frac{ρ_{i}}{ρ} (\nabla_{β} σ_{β α}^{(c)} - \nabla_{α} p + \nabla_{β} (Q_{α β γ ε} \nabla_{γ} v_{ε}))$ in lieu of the divergence of a partial stress term $\nabla_{β} T_{β α}^{(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),¹⁷ 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:¹⁷ “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

\begin{matrix} 0 = - s \nabla_{α} T - \sum_{i = 1}^{n} (α_{i 0} (v_{α}^{(i)} - \frac{1}{n} (\sum_{k = 1}^{n} v_{α}^{(k)}))) - α_{00} \frac{q_{α}}{s T} = - s \nabla_{α} T - \sum_{i = 1}^{n} (α_{i 0} v_{α}^{(i)}) - α_{00} \frac{q_{α}}{s T}, \\ 0 = \begin{matrix} \begin{matrix} - \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} (\frac{μ_{i}}{M_{i}} - \frac{μ_{j}}{M_{j}}) - α_{i 0} \frac{q_{α}}{s T} - \sum_{j = 1}^{n} (α_{i j} + α_{j i}) (v_{α}^{(i)} - v_{α}^{(j)}), \end{matrix} & i = 1, 2, \dots, n, \end{matrix} \end{matrix}

(37)

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, ${\underline{v}}^{(i)}$ by the relative ones ${\overset{⌢}{\underline{v}}}^{(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

q_{α} = - λ \nabla_{α} T + R_{G} T \sum_{\begin{matrix} i = 1 \\ j = 1 \end{matrix}}^{n} \frac{x_{j} D_{i}^{T}}{M_{i} D_{i j}} (v_{α}^{(i)} - v_{α}^{(j)}) = - λ \nabla_{α} T + R_{G} T \sum_{i = 1}^{n} (\sum_{j = 1}^{n} \frac{x_{i} x_{j}}{D_{i j}} (\frac{D_{i}^{T}}{M_{i} x_{i}} - \frac{D_{j}^{T}}{M_{j} x_{j}})) v_{α}^{(i)},

(38)

where $λ$ is the thermal conductivity, and the $D_{i j} and D_{i}^{T}$ 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 $α_{i 0}, i = 0, 1, \dots, n$ ⁠:

α_{00} = \frac{s^{2} T}{λ}, α_{i 0} = - \frac{s R_{G} T}{λ} (\sum_{j = 1}^{n} \frac{x_{i} x_{j}}{D_{i j}} (\frac{D_{i}^{T}}{M_{i} x_{i}} - \frac{D_{j}^{T}}{M_{j} x_{j}})) .

(39)

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

\sum_{j = 1}^{n} \frac{x_{i} x_{j}}{D_{i j}} (v_{α}^{(j)} - v_{α}^{(i)}) = d_{α}^{(i)} + \frac{1}{c T} \nabla_{α} T (\sum_{j = 1}^{n} \frac{x_{i} x_{j}}{D_{i j}} (\frac{D_{i}^{T}}{M_{i} x_{i}} - \frac{D_{j}^{T}}{M_{j} x_{j}})),

(40)

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

{\underline{d}}^{(j)} = \underline{\nabla} x_{j} + (x_{j} - \frac{ρ_{j}}{ρ}) \underline{\nabla} (ln p) .

(41)

Now, for an ideal multicomponent gas, the chemical potential of component j, $μ_{j}$ ⁠, is given as [see, for example, Eq. (6.1.8) of Ref. 61]

μ_{j} = μ_{0 j} (T) + R_{G} T (ln p + ln x_{j}) ideal gas .

(42)

Therefore, with the use of Eq. (42), the gradient of the chemical potential, $\underline{\nabla} μ_{j}$ at constant temperature, with temperature (or, equivalently, entropy) being an independent variable of the problem, is given as

\underline{\nabla} μ_{j} = R_{G} T (\underline{\nabla} ln (x_{j}) + \underline{\nabla} ln p) = \frac{R_{G} T}{x_{j}} (\underline{\nabla} x_{j} + x_{j} \underline{\nabla} ln p) ideal gas, constant temperature,

(43)

leading to the following equality:

\begin{matrix} \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} (\frac{μ_{i}}{M_{i}} - \frac{μ_{j}}{M_{j}}) = x_{i} c \nabla_{α} μ_{i} - \frac{ρ_{i}}{ρ} c \sum_{j = 1}^{n} x_{j} \nabla_{α} μ_{j} \\ = c R_{G} T (\underline{\nabla} x_{j} + x_{j} \underline{\nabla} ln p - \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} (\underline{\nabla} x_{j} + x_{j} \underline{\nabla} ln p)) \\ = c R_{G} T (\underline{\nabla} x_{j} + x_{j} \underline{\nabla} ln p - \frac{ρ_{i}}{ρ} \underline{\nabla} ln p) = c R_{G} T {\underline{d}}^{(j)}, ideal gas, constant temparature . \end{matrix}

(44)

Using Eq. (44), and the known value for the parameter $α_{i 0}$ as provided by Eq. (39), Eq. (40) becomes

\frac{1}{c R_{G} T} \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} (\frac{μ_{i}}{M_{i}} - \frac{μ_{j}}{M_{j}}) = \sum_{j = 1}^{n} \frac{x_{i} x_{j}}{D_{i j}} (v_{α}^{(j)} - v_{α}^{(i)}) + \frac{1}{c T} \nabla_{α} T \frac{α_{i 0} λ}{s R_{G} T} .

(45)

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

\frac{1}{c R_{G} T} \frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} (\frac{μ_{i}}{M_{i}} - \frac{μ_{j}}{M_{j}}) = \sum_{j = 1}^{n} \frac{x_{i} x_{j}}{D_{i j}} (v_{α}^{(j)} - v_{α}^{(i)}) - \frac{α_{i 0} λ}{c R_{G} s^{2} T^{2}} (\sum_{i = 1}^{n} (α_{i 0} v_{α}^{(i)}) + α_{00} \frac{q_{α}}{s T}),

(46)

or, using Eq. (39) to eliminate $λ$ and after multiplying both sides by $c R_{G} T$ we find

\frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} (\frac{μ_{i}}{M_{i}} - \frac{μ_{j}}{M_{j}}) = c R_{G} T \sum_{j = 1}^{n} \frac{x_{i} x_{j}}{D_{i j}} (v_{α}^{(j)} - v_{α}^{(i)}) - \frac{α_{i 0}}{α_{00}} (\sum_{j = 1}^{n} (α_{j 0} v_{α}^{(j)})) - α_{i 0} \frac{q_{α}}{s T} .

(47)

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

\frac{ρ_{i}}{ρ} \sum_{j = 1}^{n} ρ_{j} \nabla_{α} (\frac{μ_{i}}{M_{i}} - \frac{μ_{j}}{M_{j}}) = c R_{G} T \sum_{j = 1}^{n} \frac{x_{i} x_{j}}{D_{i j}} (v_{α}^{(j)} - v_{α}^{(i)}) - \frac{α_{i 0}}{α_{00}} (\sum_{j = 1}^{n} (α_{j 0} (v_{α}^{(j)} - v_{α}^{(i)}))) - α_{i 0} \frac{q_{α}}{s T} .

(48)

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 $α_{i j}$ for dilute gases:¹²

α_{i j} + α_{j i} = c R_{G} T \frac{x_{i} x_{j}}{D_{i j}} - \frac{α_{i 0} α_{j 0}}{α_{00}} dilute gases .

(49)

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¹² 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), $α_{10} = - α_{20}$ (as mentioned before), and the only nonzero mass flux coefficient following Eq. (49) is $α_{12}$ ⁠, which is always non-negative:

2 α_{12} = 2 {\hat{α}}_{12} + \frac{α_{10}^{2}}{α_{00}},

(50)

where

2 {\hat{α}}_{12} \equiv c R_{G} T \frac{x_{1} x_{2}}{D_{12}} \geq 0.

(51)

For dilute binary gas, the significance of the reduced coefficient ${\hat{α}}_{i j}$ 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 $α_{12}, α_{00}, μ, κ^{'}$ being non-negative

2 α_{12} α_{00} - {(α_{10})}^{2} \geq 0 \Leftrightarrow α_{00} (2 α_{12} - \frac{{(α_{10})}^{2}}{α_{00}}) \geq 0 \Leftrightarrow 2 α_{00} {\hat{α}}_{12} \geq 0 dilute gases,

(52)

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 ${\hat{α}}_{12}$ ⁠. 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, $N_{2} - H_{2}$ ⁠, 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 $λ = 2.0 \times 10^{- 4} cal / (cm s ° K)$ ⁠, whereas in Table 8.4-12 a binary diffusivity at the same temperature is reported as $D_{12} = 0.67 {cm}^{2} s^{- 1}$ ⁠. 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 $k_{T} \approx 0.068$ ⁠. From the definition of the thermal diffusion ratio for a binary gas mixture [Eq. (8.1-13) of Ref. 12]

{\underline{v}}^{(1)} - {\underline{v}}^{(2)} = - \frac{1}{x_{1} x_{2}} D_{12} (\underline{\nabla} x_{1} + k_{T} \frac{1}{T} \underline{\nabla} T) .

(53)

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:

k_{T} = \frac{1}{c} (\frac{x_{1} x_{2}}{D_{12}} (\frac{D_{1}^{T}}{M_{1} x_{1}} - \frac{D_{2}^{T}}{M_{2} x_{2}})) .

(54)

This results in the following relationship for the $α_{10}$ coefficient:

α_{10} = - \frac{s R_{G} T}{λ} (\frac{x_{1} x_{2}}{D_{12}} (\frac{D_{1}^{T}}{M_{1} x_{1}} - \frac{D_{2}^{T}}{M_{2} x_{2}})) = - \frac{s R_{G} T}{λ} c k_{T} .

(55)

Finally, the total molar concentration for an ideal gas at the same temperature is $1 / 22 400 mol {cm}^{- 3}$ ⁠. With all this information, it is straightforward to evaluate the relative correction to the reduced mass transfer coefficient ${\hat{α}}_{12}$

\frac{α_{10}^{2}}{2 {\hat{α}}_{12} α_{00}} = \frac{{(\frac{s R_{G} T}{λ} c k_{T})}^{2}}{c R_{G} T \frac{x_{1} x_{2}}{D_{12}} \frac{s^{2} T}{λ}} = \frac{R_{G} c k_{T}^{2} D_{12}}{x_{1} x_{2} λ} = \frac{2 \times {(0.068)}^{2} 0.67}{0.65 \times 0.35 \times 0.0002 \times 22400} = 0.0061 .

(56)

We therefore see that at least in this case this correction is small and can be neglected, thereby approximating $α_{12} \approx {\hat{α}}_{12}$ ⁠, 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

\begin{matrix} x_{A} ≪ 1; x_{B} = 1 - x_{A} \approx 1; ρ \approx constant \approx ρ_{B} \approx c M_{B}; w_{A} = M_{A} c_{A} = M_{A} c x_{A} \\ μ_{A} = g_{A} (T, P) + R_{G} T ln x_{A}; μ_{B} \approx g_{B} (T, P); \underline{v} \approx {\underline{v}}_{B} \Rightarrow {\overset{⌢}{v}}_{- B} \approx \underline{0} \end{matrix}

(57)

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, j = B$ ⁠:

\begin{matrix} \frac{ρ_{A}}{ρ} (ρ_{B} \nabla_{α} (\frac{μ_{A}}{M_{A}} - \frac{μ_{B}}{M_{B}})) = R_{G} T \frac{ρ_{B}}{ρ_{B}} \frac{M_{A} x_{A} c}{M_{A} x_{A}} \nabla_{α} x_{A} = R_{G} T c \nabla_{α} x_{A} = \frac{R_{G} T ρ}{M_{A}} \nabla_{α} w_{A}, \\ (c R_{G} T \frac{x_{A} x_{B}}{D_{AB}}) (v_{B α} - v_{A α}) = c R_{G} T \frac{x_{A}}{D_{AB}} (- \frac{j_{A α}}{ρ_{A}}) = c R_{G} T \frac{1}{M_{A} c D_{AB}} (- j_{A α}) = - \frac{R_{G} T}{M_{A} D_{AB}} j_{A α}, \end{matrix}

(58)

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

j_{A α} = - ρ D_{AB} \nabla_{α} w_{A},

(59)

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

\begin{matrix} \overset{Δ}{(ρ_{A} {\overset{⌢}{v}}_{A α})} = - ρ_{A} \nabla_{α} (- \frac{1}{2} {({\overset{⌢}{v}}_{A})}^{2} + \frac{μ_{A}}{M_{i}} - \frac{μ_{B}}{M_{j}}) - 2 α_{AB} ({\overset{⌢}{v}}_{A α}) \\ = - \frac{R_{G} T}{M_{A}} \nabla_{α} ρ_{A} + ρ_{A} \nabla_{α} (\frac{1}{2} {({\overset{⌢}{v}}_{A})}^{2}) - \frac{R_{G} T}{M_{A}} \frac{ρ_{A}}{D_{AB}} ({\overset{⌢}{v}}_{A α}) . \end{matrix}

(60)

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

\begin{matrix} α_{10} (\frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(1)}} - \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(2)}}) + α_{00} \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(0)}} = - s \underline{\nabla} \frac{δ H}{δ s}, \\ (α_{12} + α_{21}) (\frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(1)}} - \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(2)}}) + α_{10} \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(0)}} = \begin{matrix} \begin{matrix} - \frac{ρ_{1} ρ_{2}}{ρ} \underline{\nabla} (\frac{δ H}{δ ρ_{1}} - \frac{δ H}{δ ρ_{2}}), \end{matrix} \end{matrix} \end{matrix}

(61)

leading to the following relations:

[\begin{matrix} (\frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(1)}} - \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(2)}}) \\ \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(0)}} \end{matrix}] = - \frac{1}{2 α_{00} {\hat{α}}_{12}} [\begin{matrix} α_{00} & - α_{10} \\ - α_{10} & 2 α_{12} \end{matrix}] [\begin{matrix} \frac{ρ_{1} ρ_{2}}{ρ} \underline{\nabla} (\frac{δ H}{δ ρ_{1}} - \frac{δ H}{δ ρ_{2}}) \\ s \underline{\nabla} \frac{δ H}{δ s} \end{matrix}] .

(62)

If we now use the suggested replacements by Eq. (55), Eqs. (51) (with the approximation $α_{12} \approx {\hat{α}}_{12}$ ⁠) and (39) we get

[\begin{matrix} (\frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(1)}} - \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(2)}}) \\ \frac{δ H}{δ {\underline{\overset{⌢}{u}}}^{(0)}} \end{matrix}] = - [\begin{matrix} \frac{D_{12}}{c R_{G} T x_{1} x_{2}} & - \frac{1}{sTc} (\frac{D_{1}^{T}}{M_{1} x_{1}} - \frac{D_{2}^{T}}{M_{2} x_{2}}) \\ - \frac{1}{sTc} (\frac{D_{1}^{T}}{M_{1} x_{1}} - \frac{D_{2}^{T}}{M_{2} x_{2}}) & \frac{λ}{s^{2} T} \end{matrix}] [\begin{matrix} \frac{ρ_{1} ρ_{2}}{ρ} \underline{\nabla} (\frac{δ H}{δ ρ_{1}} - \frac{δ H}{δ ρ_{2}}) \\ s \underline{\nabla} \frac{δ H}{δ s} \end{matrix}] .

(63)

This results in the following symmetric dissipation bracket:

{[F, H]}_{s p} = - \int_{Ω} {[\begin{matrix} \underline{\nabla} (\frac{δ F}{δ ρ_{1}} - \frac{δ F}{δ ρ_{2}}) \\ \underline{\nabla} \frac{δ F}{δ s} \end{matrix}]}^{T} \cdot [\begin{matrix} \frac{D_{12} c_{1} c_{2}}{c R_{G} T} & - \frac{1}{T} (ρ_{2} \frac{D_{1}^{T}}{ρ} - ρ_{1} \frac{D_{2}^{T}}{ρ}) \\ - \frac{1}{T} (ρ_{2} \frac{D_{1}^{T}}{ρ} - ρ_{1} \frac{D_{2}^{T}}{ρ}) & \frac{λ}{T} \end{matrix}] \cdot [\begin{matrix} \underline{\nabla} (\frac{δ H}{δ ρ_{1}} - \frac{δ H}{δ ρ_{2}}) \\ \underline{\nabla} \frac{δ H}{δ s} \end{matrix}] d^{3} x .

(64)

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., $D_{12}, λ, D_{1}^{T}, D_{2}^{T}$ ⁠) with explicit expressions also provided by kinetic theory.¹² Even if these correlations have been derived for binary mixtures before⁶⁵ 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 $α_{12}$ and ${\hat{α}}_{12}$ ⁠). 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

\frac{\partial ρ_{A}}{\partial t} = - \nabla_{β} (v_{A β} ρ_{A}),

(65)

and (2) the corresponding species momentum equation, Eq. (60) for a stagnant medium, for which $\underline{v} \approx {\underline{v}}_{B} = \underline{0} \Rightarrow {\underline{\overset{⌢}{v}}}_{A} = {\underline{v}}_{A}$ and further, by neglecting higher order contributions, yields the following linearized equation:

\frac{\partial (ρ_{A} v_{A β})}{\partial t} = - \frac{R_{G} T ρ}{M_{A}} \nabla_{β} w_{A} - \frac{R_{G} T}{M_{A} D_{AB}} (ρ_{A} v_{A β}) .

(66)

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

ρ \frac{\partial^{2} (w_{A})}{\partial t^{2}} = \frac{R_{G} T ρ}{M_{A}} \nabla^{2} w_{A} + \frac{R_{G} T}{M_{A} D_{AB}} \nabla_{β} (ρ_{A} v_{A β}) \Leftrightarrow \frac{\partial^{2} (w_{A})}{\partial t^{2}} = \frac{R_{G} T}{M_{A}} \nabla^{2} w_{A} - \frac{R_{G} T}{M_{A} D_{AB}} \frac{\partial (w_{A})}{\partial t},

(67)

with characteristic speed

c = \pm \sqrt{\frac{R_{G} T}{M_{A}}},

(68)

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, $M_{A} = 18 kg / kmol$ ⁠) and for room temperature, $T = 298 K$

c = \pm \sqrt{\frac{R_{G} T}{M_{A}}} = \pm \sqrt{\frac{8314.51 \times 298}{18}} m / s = 371 m / s,

(69)

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 $v_{rms} = \sqrt{\frac{(3 R_{G} T)}{M_{A}}}$ ⁠, enables to rewrite Eq. (67) as a hyperbolic equation

\frac{1}{c^{2}} \frac{\partial^{2} (w_{A})}{\partial t^{2}} = \nabla^{2} w_{A} - \frac{1}{D_{AB}} \frac{\partial (w_{A})}{\partial t},

(70)

which reduces to the usual, linear parabolic diffusion equation when $c^{2} \to \infty$ ⁠. 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⁶⁶ and later in heat transfer by Cattaneo.⁶⁷ 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.⁶⁹ 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

\frac{\partial w_{A}}{\partial t} + \frac{D_{AB}}{c^{2}} \frac{\partial^{2} w_{A}}{\partial t^{2}} = D_{AB} \nabla^{2} w_{A} \Leftrightarrow \frac{\partial w_{A}}{\partial t} + τ \frac{\partial^{2} w_{A}}{\partial t^{2}} = D_{AB} \nabla^{2} w_{A} .

(71)

Here, $τ$ is the relaxation time scale for propagation defined as $τ = \frac{D_{AB}}{c^{2}}$ 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⁶⁷ to account for inertia arising from finite speed of propagation in heat transfer.

This problem was further explored by Sieniutycz and Berry⁷⁰ 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⁷¹ 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 \times 10^{- 8}$ to $3 \times 10^{- 2}$ 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 $2 L$ ⁠, where the solute enters through both boundaries (applied as a step function) at $t = 0$ with weight fraction $w_{A}^{(1)}$ ⁠, with the solution initially being at solute weight fraction of $w_{A}^{(0)}$ ⁠. Let us further consider as spatial domain of interest the interval $[- 2 L, 0]$ ⁠. For this system, the governing equation can be expressed by non-dimensionalizing Eq. (71) using the weight fraction at boundary $w_{A}^{(1)}$ ⁠, the initial weight fraction $w_{A}^{(0)}$ ⁠, the characteristic time $τ$ ⁠, and the corresponding diffusion length $\sqrt{D_{AB} τ}$ as scales to construct appropriate dimensionless variables (denoted here with * superscripts), such that $x^{*} = \frac{x}{\sqrt{D_{AB} τ}}$ ⁠, $L^{*} = \frac{L}{\sqrt{D_{AB} τ}}$ ⁠, $w_{A}^{*} = \frac{w_{A} - w_{A}^{(1)}}{w_{A}^{(0)} - w_{A}^{(1)}}$ and $t^{*} = t / τ$ ⁠.

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^{*} \leq O (1)$ ⁠. Sandler and Dahler⁵⁹ argue that the differences are only significant at length scales of $100 Å$ in liquids or $10^{5} Å$ in gases.

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

\frac{\partial w_{A}^{*}}{\partial t^{*}} + \frac{\partial^{2} w_{A}^{*}}{\partial t^{*}^{2}} = \frac{\partial^{2} w_{A}^{*}}{\partial x^{*}^{2}}, and w_{A}^{*} (x^{*}, t^{*} = 0) = 1; w_{A}^{*} (x^{*} = - 2 L^{*}, t^{*}) = 0; w_{A}^{*} (x^{*} = 0, t^{*}) = 0.

(72)

Exploring the symmetry over the centerline, this problem can be solved in half domain $(- L^{*}, 0)$ and mirrored about $x^{*} = - L^{*}$ ⁠. We can therefore replace the left Dirichlet boundary condition at $x^{*} = - 2 L^{*}$ with a Neumann boundary condition at $x^{*} = - L^{*}$ such that $\frac{\partial w_{A}^{*}}{\partial x^{*}} (x^{*} = - L^{*}, t^{*}) = 0$ ⁠. One additional boundary condition is required to obtain the coefficients, since the equation in second order in time, for which we propose $\frac{\partial w_{A}^{*}}{\partial t^{*}} (x^{*}, t^{*} = 0) = 0$ ⁠. Solving the above equation using Fourier eigenfunctions such that $w_{A}^{*} = \sum_{n = 0}^{\infty} Θ_{n} (t^{*}) Φ_{n} (x^{*})$ ⁠, the eigenfunctions and eigenvalues, $Φ_{n} (x^{*}), λ_{n}$ ⁠, respectively, based on this combination of homogeneous Neumann and Dirichlet boundary conditions come out to be

\begin{matrix} Φ_{n} (x^{*}) = \sqrt{\frac{2}{L^{*}}} cos (n + \frac{1}{2}) \frac{π}{L^{*}} (x^{*} + L^{*}) = \sqrt{\frac{2}{L^{*}}} cos λ_{n} (x^{*} + L^{*}), \\ λ_{n} = (n + \frac{1}{2}) \frac{π}{L^{*}}, \end{matrix}

(73)

where $n = 0, 1, 2, \dots$ ⁠. Using a linear projection of variables along the eigenfunctions, defined as $Θ_{n} (t^{*}) = \int_{0}^{L^{*}} Φ_{n} (x^{*}) w_{A}^{*} d x$ ⁠, the transformed equation and boundary conditions become

\frac{\partial Θ_{n}}{\partial t^{*}} + \frac{\partial^{2} Θ_{n}}{\partial t^{*}^{2}} + λ_{n}^{2} τ^{*} Θ_{n} = 0, with Θ_{n} (x^{*}, t^{*} = 0) = \sqrt{\frac{2}{L^{*}}} \frac{{(- 1)}^{n}}{λ_{n}} and \frac{\partial Θ}{\partial t^{*}} (x^{*}, t^{*} = 0) = 0.

(74)

The solution for this transformed equation is

Θ_{n} (t^{*}) = \sqrt{\frac{2}{L^{*}}} e^{- \frac{t^{*}}{2}} {(\frac{{(- 1)}^{n}}{2 λ_{n}} - \frac{{(- 1)}^{n}}{4 λ_{n} β}) e^{β t^{*}} + (\frac{{(- 1)}^{n}}{2 λ_{n}} + \frac{{(- 1)}^{n}}{4 λ_{n} β}) e^{- β t^{*}}},

(75)

where $β = \frac{\sqrt{1 - 4 λ_{n}^{2} τ^{*}}}{2}$ ⁠. Combining Eqs. (73) and (75), we can construct the solution for the dimensionless variable $w_{A}^{*}$ as follows:

w_{A}^{*} = \frac{2}{L^{*}} \sum_{n = 0}^{\infty} e^{- \frac{t^{*}}{2}} {(\frac{{(- 1)}^{n}}{2 λ_{n}} - \frac{{(- 1)}^{n}}{4 λ_{n} β}) e^{- β t^{*}} + (\frac{{(- 1)}^{n}}{2 λ_{n}} + \frac{{(- 1)}^{n}}{4 λ_{n} β}) e^{β t^{*}}} cos λ_{n} (x^{*} + L^{*}),

(76)

which can be expressed in the dimensional space as

w_{A} = w_{A}^{(1)} - \frac{2 \sqrt{D_{AB} τ}}{L} (w_{A}^{(1)} - w_{A}^{(0)}) \sum_{n = 0}^{\infty} e^{- \frac{t}{2 τ}} {(\frac{{(- 1)}^{n}}{2 λ_{n}} - \frac{{(- 1)}^{n}}{4 λ_{n} β}) e^{- β \frac{t}{τ}} + (\frac{{(- 1)}^{n}}{2 λ_{n}} + \frac{{(- 1)}^{n}}{4 λ_{n} β}) e^{β \frac{t}{τ}}} cos (n + \frac{1}{2}) π \frac{x + L}{\sqrt{D_{AB} τ}} .

(77)

Figure 1 shows the above solution plotted for an initial solute weight fraction $w_{A}^{(0)} = 0$ and weight fraction at the boundary $w_{A}^{(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.

FIG. 1.

View large Download slide

Series solution for transient diffusion problem comparing the parabolic and hyperbolic equations. The problem was solved in a one-dimensional slab of length of length $2 L$ ⁠, with $L = 10 \sqrt{D_{AB} τ}$ ⁠. At short times, a sharp discontinuity is seen in the species front which is consistent with the finite speed of propagation. On the other hand, the parabolic diffusion equation or Fick's law overpredicts the concentration of species at short times, but as the solution approaches equilibrium, the two solutions become virtually indistinguishable.

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 (⁠ $- \infty \leq x^{*} \leq \infty$ ⁠) of an initial delta function, i.e., subject to the initial conditions,

w_{A} (x^{*}, 0) = δ (x^{*}); \frac{\partial w_{A}}{\partial t^{*}} (x^{*}, 0) = 0.

(78)

The solution to this problem has been developed by Gomez⁷² based on modified Bessel functions of the first kind. The solution was found to be

w_{A} (x^{*}, t^{*}) \equiv w_{A}^{G} (x^{*}, t^{*}) = {\begin{array}{l} \frac{1}{2} e^{- \frac{t^{*}}{2}} [δ (| x^{*} | - t^{*}) + \frac{1}{2} I_{0} (z) + \frac{1}{4} t^{*} \frac{I_{1} (z)}{z}], & | x^{*} | \leq t^{*} \\ 0, & | x^{*} | > t^{*} \end{array} z \equiv \frac{1}{2} \sqrt{t^{*}^{2} - x^{*}^{2}},

(79)

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

I_{ν} (x) \equiv i^{- ν} J_{ν} (i x) = \sum_{k = 0}^{\infty} \frac{x^{2 k + ν}}{2^{2 k + ν} (k!) ((k + v)!)}, v = 0, 1, 2, \dots .

(80)

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

w_{A} (x^{*}, 0) = f (x^{*}); \frac{\partial w_{A}}{\partial t^{*}} (x^{*}, 0) = 0,

(81)

can therefore be obtained as a convolution

w_{A} (x^{*}, t^{*}) = \int_{- \infty}^{\infty} f (x') w_{A}^{G} (x^{*} - x', t^{*}) d x', t^{*} \geq 0.

(82)

For the semi-infinite domain (⁠ $x^{*} \leq 0$ ⁠) and the Dirichlet boundary condition applied on the right boundary [ $w_{A} (x^{*} = 0, t^{*}) = 1$ ], the corresponding f function is

f (x') = {\begin{matrix} 2, & x' \geq 0, \\ 0, & x' < 0. \end{matrix}

(83)

Furthermore, because for all values of $x^{*} \leq - t^{*}$ the solution is equal to 0, Eq. (82) can be simplified to

w_{A} (x^{*}, t^{*}) = \int_{0}^{max (0, x^{*} + t^{*})} f (x') w_{A}^{G} (x^{*} - x', t^{*}) d x', x^{*} \leq 0; t^{*} \geq 0.

(84)

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

\begin{matrix} w_{A} (x^{*}, t^{*}) = 2 \int_{- t^{*}}^{x^{*}} \frac{1}{2} e^{- \frac{t^{*}}{2}} [δ (| \hat{x} | - t^{*}) + \frac{1}{2} I_{0} (\hat{z}) + \frac{1}{4} t^{*} \frac{I_{1} (\hat{z})}{\hat{z}}] d \hat{x}, - t^{*} \leq \hat{x} \leq 0; \hat{z} \equiv \frac{1}{2} \sqrt{t^{* 2} - {\hat{x}}^{2}} \\ = e^{- \frac{t^{*}}{2}} + \frac{1}{2} e^{- \frac{t^{*}}{2}} [\int_{- t^{*}}^{x^{*}} I_{0} (\hat{z}) d \hat{x}] + \frac{1}{4} t^{*} e^{- \frac{t^{*}}{2}} [\int_{- t^{*}}^{x^{*}} \frac{I_{1} (\hat{z})}{\hat{z}} d \hat{x}], - t^{*} \leq \hat{x} \leq 0; \hat{z} \equiv \frac{1}{2} \sqrt{t^{* 2} - {\hat{x}}^{2}} . \end{matrix}

(85)

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^{*} = - t^{*}$ ⁠. 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 (⁠ $Δ t^{*}$ ⁠) to be smaller or equal than the spatial steps (⁠ $Δ x^{*}$ ⁠). For simplicity and for illustration purposes, we propose here a first order finite difference scheme with the time steps (⁠ $Δ t^{*}$ ⁠) equal to the spatial steps (⁠ $Δ 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:

\begin{matrix} D_{+} \equiv \frac{\partial}{\partial t^{*}} - \frac{\partial}{\partial x^{*}}, \\ D_{-} \equiv \frac{\partial}{\partial t^{*}} + \frac{\partial}{\partial x^{*}}, \end{matrix}

(86)

where $D_{+}$ is computed along the $x^{*} = t^{*}$ (positive) characteristic, and $D_{-}$ is computed along the $x^{*} = - t^{*}$ (negative) characteristic. Equation (72) can then be written in terms of these operators as

D_{+} D_{-} w_{A} + \frac{1}{2} (D_{+} + D_{-}) w_{A} = 0.

(87)

FIG. 2.

View large Download slide

Discretization scheme used for the finite difference solution of the transient diffusion equation. The filled circles indicate staggered grid points located at half steps between two grid points (indicated by open circles). The spatial coordinates use index j, and the time steps use index i.

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

\begin{matrix} D_{+} D_{-} w_{A} = (\frac{\frac{w_{A}^{(i, j)} - w_{A}^{(i - \frac{1}{2}, j - \frac{1}{2})}}{\frac{Δ t^{*}}{2}} - \frac{w_{A}^{(i - \frac{1}{2}, j + \frac{1}{2})} - w_{A}^{(i - 1, j)}}{\frac{Δ t^{*}}{2}}}{\frac{Δ t^{*}}{2}}), \\ \frac{1}{2} (D_{+} + D_{-}) w_{A} = \frac{w_{A}^{(i, j)} - w_{A}^{(i - 1, j)}}{Δ t^{*}} . \end{matrix}

(88)

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:

D^{-} Δ w + \frac{1}{2} Δ w = 0,

(89)

which readily admits the analytical solution $Δ w = exp (- \frac{t}{2})$ ⁠. 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^{*} \leq 0$ ⁠) 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 $τ$ ⁠.

FIG. 3.

View large Download slide

Comparison of analytical solution using the modified Bessel function of the first kind and finite differences. The solutions are presented with three different spatial resolutions. The graph at the bottom indicates deviation of the finite difference solutions from the analytic solution.

FIG. 4.

$FIG. 4. Comparison of finite differences with built-in MATLAB parabolic-elliptic PDE solver, PDEPE. Both parabolic and hyperbolic transient diffusion equations are solved in a one-dimensional slab of length 2L, with L=10DABτ. Dirichlet boundary condition wA=1 is applied on both boundaries, with initial weight fraction of species A, wA(0)=0 everywhere in the domain.$

View large Download slide

Comparison of finite differences with built-in MATLAB parabolic-elliptic PDE solver, PDEPE. Both parabolic and hyperbolic transient diffusion equations are solved in a one-dimensional slab of length $2 L$ ⁠, with $L = 10 \sqrt{D_{A B} τ}$ ⁠. Dirichlet boundary condition $w_{A} = 1$ is applied on both boundaries, with initial weight fraction of species $A$ ⁠, $w_{A}^{(0)} = 0$ everywhere in the domain.

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.³⁷ 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^{*} ≫ 1$ ⁠, as it allows for the use of efficient hyperbolic numerical solvers. Of course, when $L^{*} ≫ 1$ 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,²¹ 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,¹⁵ 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

$a_{j}$: $j th$ species/component in the system
$c$: molar concentration
$c_{s}$: characteristic speed
${\underline{\underline{C}}}^{(k)}$: conformation tensor corresponding to component $k$
$d V$ ⁠, $d^{3} x$: differential volume
$D_{i j}$: molecular diffusion coefficient between $i th - j th$ species
$D_{i}^{T}$: Thermal diffusion coefficient of $i th$ species
$h$: Hamiltonian functional density
$H$: Hamiltonian functional
$\underline{\underline{I}}$: identity tensor
$\underline{j}$: mass-based diffusive flux of $j th$ species
$k_{T}$: thermal diffusion ratio
$L$: length
$M_{j}$: molar mass/molecular weight
$p$: thermodynamic pressure
$\underline{q}$: heat flux
$\underline{\underline{Q}}$: viscous dissipation tensor
$R_{G}$: universal gas constant
$s$: total entropy density
$t$: time
$T$: temperature of the system
$\underline{u}$: total momentum density
${\underline{u}}^{(j)}$: $j th$ species momentum density, $j \neq 0$
${\underline{\overset{⌢}{u}}}^{(j)}$: $j th$ species reduced momentum density, $j \neq 0$
${\underline{\overset{⌢}{u}}}^{(0)}$: heat momentum density
$\underline{v}$: velocity field of the system or reference frame
${\underline{v}}^{(j)}$: $j th$ species velocity, $j \neq 0$
${\underline{\overset{⌢}{v}}}^{(j)}$: $j th$ species reduced velocity, $j \neq 0$
${\underline{\overset{⌢}{v}}}^{(0)}$: heat velocity
$w_{A}$: weight fraction of species A
$x$: spatial coordinate
$x_{j}$: molar fraction of $j th$ species

Greek symbols

$α_{i j}$: generalized coefficient of mass transfer between components $i$ and $j$ ⁠, $i, j \neq 0$
${\hat{α}}_{i j}$: reduced mass transfer coefficient, $i, j \neq 0$
$α_{i 0}$: Soret–Dufour coefficient for component $i$ ⁠, $i \neq 0$
$α_{00}$: generalized heat conductivity coefficient
$ε$: energy density functional
$κ'$: bulk viscosity
$λ$: thermal conductivity
$μ$: shear viscosity
$μ_{j}$: $j th$ species chemical potential
$μ_{0 j}$: reference chemical potential of $j th$ species
$Ξ$: primitive dissipation functional
$ρ$: total system density
$ρ_{j}$: $j th$ species density
$ρ_{0}$: effective heat mass density
$\underline{\underline{σ}}$: stress tensor
$σ_{s}$: volumetric rate of entropy production
$τ$: characteristic time in hyperbolic diffusion
$Ω$: domain of interest in $\underline{x}$ space

Brackets/scripts

$[\cdot, \cdot]$: dissipation bracket
${[\cdot, \cdot]}_{p}$: primitive dissipation bracket
${[\cdot, \cdot]}_{e c}$: entropy correction term
${[\cdot, \cdot]}_{a}$: anti-symmetric part of dissipation bracket
${[\cdot, \cdot]}_{s p}$: symmetric part of dissipation bracket
${\cdot, \cdot}$: Poisson bracket
$⟨ \cdot ⟩$: average/mean
$\underline{\overset{Δ}{\cdot}}, \overset{Δ}{\underline{\underline{\cdot}}}$: lower-convected Oldroyd time derivative for a tensor
$\vec{\cdot}$: quantity in n-vector space
$\vec{\vec{\cdot}}$: matrix in $n \times n$ 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 $\underline{w}$ and $\underline{\overset{⌢}{w}}$ 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:

{\underline{\overset{⌢}{u}}}^{(j)} = {\underline{u}}^{(j)} - \frac{ρ_{j}}{ρ} \sum_{k = 1}^{n} {\underline{u}}^{(k)}; \underline{u} = \sum_{k = 1}^{n} {\underline{u}}^{(k)}; subject to the constraint \sum_{k = 1}^{n} {\underline{\overset{⌢}{u}}}^{(k)} = \underline{0} .

(A1)

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

{\underline{u}}^{(j)} = {\underline{\overset{⌢}{u}}}^{(j)} + \frac{ρ_{j}}{ρ} \underline{u} .

(A2)

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:

\frac{δ G}{δ {\underline{u}}^{(j)}} = \frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(j)}} - \frac{1}{ρ} \sum_{k = 1}^{n} (ρ_{k} \frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(k)}}) + \frac{δ G}{δ \underline{u}}

(A3)

and

\frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(j)}} = \frac{δ G}{δ {\underline{u}}^{(j)}} - \frac{1}{n} \sum_{k = 1}^{n} \frac{δ G}{δ {\underline{u}}^{(k)}}; \frac{δ G}{δ \underline{u}} = \frac{1}{ρ} \sum_{k = 1}^{n} (ρ_{k} \frac{δ G}{δ {\underline{u}}^{(k)}}) .

(A4)

Note the correction term to the first relationship in Eq. (A4) due to the presence of the constraint, Eq. (7).

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

Considering the variables $\underline{\overset{⌢}{w}} \equiv {[ρ_{1}, ρ_{2}, \dots, ρ_{n}, \underline{u}, {\underline{\overset{⌢}{u}}}^{(1)}, {\underline{\overset{⌢}{u}}}^{(2)}, \dots, {\underline{\overset{⌢}{u}}}^{(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

\begin{matrix} {[F, H]}_{s p} = - \int_{Ω} \sum_{i, j = 1}^{n} a_{i j} \frac{δ F}{δ {\overset{⌢}{u}}_{α}^{(i)}} \frac{δ H}{δ {\overset{⌢}{u}}_{α}^{(j)}} d^{3} x \\ - \int_{Ω} Q_{α β γ ε} \nabla_{α} \frac{δ F}{δ u_{β}} \nabla_{γ} \frac{δ H}{δ u_{ε}} d^{3} x \\ - \int_{Ω} \frac{κ}{T} (\nabla_{α} \frac{δ F}{δ s}) (\nabla_{α} \frac{δ H}{δ s}) d^{3} x, \end{matrix}

(B1)

where $Q_{α β γ ε} \equiv μ (δ_{α γ} δ_{β ε} + δ_{α ε} δ_{β γ}) + κ^{'} δ_{α β} δ_{γ ε}$ with $μ, κ^{'}$ are the shear and bulk viscosity⁴¹ and $a_{i j}, κ$ 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 $a_{i j}, κ$ need to be substituted by symmetric second order tensors, while $Q_{α β γ ε}$ 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 $\vec{\underline{\overset{⌢}{u}}} \equiv {[{\underline{\overset{⌢}{u}}}^{(1)}, {\underline{\overset{⌢}{u}}}^{(2)}, \dots, {\underline{\overset{⌢}{u}}}^{(n)}]}^{T}$ with the unit n-vector $\vec{k} \equiv \frac{1}{\sqrt{n}} {[\underset{n terms}{\underset{︸}{1, 1, \dots, 1}}]}^{T}$ where the Euclidean L2 norm is used for normalization

{\vec{k}}^{T} \cdot \vec{\underline{\overset{⌢}{u}}} = \frac{1}{\sqrt{n}} \sum_{j = 1}^{n} {\underline{\overset{⌢}{u}}}^{(j)} = \underline{0},

(B2)

i.e., the n-vector $\vec{\underline{\overset{⌢}{u}}} \equiv {[{\underline{\overset{⌢}{u}}}^{(1)}, {\underline{\overset{⌢}{u}}}^{(2)}, \dots, {\underline{\overset{⌢}{u}}}^{(n)}]}^{T}$ is restricted to be in a subspace that is perpendicular to with the unit n-vector $\vec{k} \equiv \frac{1}{\sqrt{n}} {[\underset{n terms}{\underset{︸}{1, 1, \dots, 1}}]}^{T}$ ⁠. Due to this constraint, it is well known that a similar condition also applies⁴¹ to the n-vector of the Volterra derivatives with respect to the reduced momenta of any arbitrary functional G (G = H, F, ….), $\vec{\frac{δ G}{δ \underline{\overset{⌢}{u}}}} \equiv {[\frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(1)}}, \frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(2)}}, \dots, \frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(n)}}]}^{T}$

{\vec{k}}^{T} \cdot \vec{\frac{δ G}{δ \underline{\overset{⌢}{u}}}} = \frac{1}{\sqrt{n}} \sum_{j = 1}^{n} \frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(j)}} = \underline{0},

(B3)

i.e., the n-vector $\vec{\frac{δ G}{δ \underline{\overset{⌢}{u}}}} \equiv {[\frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(1)}}, \frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(2)}}, \dots, \frac{δ G}{δ {\underline{\overset{⌢}{u}}}^{(n)}}]}^{T}$ is restricted to be in a subspace that is perpendicular to the unit n-vector $\vec{k} \equiv \frac{1}{\sqrt{n}} {[\underset{n terms}{\underset{︸}{1, 1, \dots, 1}}]}^{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 $\vec{\vec{a}}$ are independent. Instead, the only relevant information (that can lead to measurable results) belongs to the subspace of the n × n coefficient symmetric matrix $\vec{\vec{a}}$ that satisfy the corresponding to Eq. (B3) condition

{\vec{k}}^{T} \cdot \vec{\vec{a}} = \vec{\vec{a}} \cdot \vec{k} = \vec{0},

(B4)

i.e., without loss of generality we can assume that the n × n coefficient symmetric matrix $\vec{\vec{a}}$ has a zero eigenvalue along the unit n-vector $\vec{k} \equiv \frac{1}{\sqrt{n}} {[\underset{n terms}{\underset{︸}{1, 1, \dots, 1}}]}^{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

\vec{\vec{a}} = [{\vec{q}}^{(1)}, {\vec{q}}^{(2)}, \dots, {\vec{q}}^{(n - 1)}, \vec{k}] [\begin{matrix} \vec{\vec{\overset{⌢}{a}}} & \vec{\overset{⌢}{0}} \\ {\vec{\overset{⌢}{0}}}^{T} & 0 \end{matrix}] {[{\vec{q}}^{(1)}, {\vec{q}}^{(2)}, \dots, {\vec{q}}^{(n - 1)}, \vec{k}]}^{T},

(B5)

where $[{\vec{q}}^{(1)}, {\vec{q}}^{(2)}, \dots, {\vec{q}}^{(n - 1)}]$ is a set of n-1 independent orthonormal n-vectors perpendicular to the unit n-vector $\vec{k} \equiv \frac{1}{\sqrt{n}} {[\underset{n terms}{\underset{︸}{1, 1, \dots, 1}}]}^{T}$ ⁠, $\vec{\overset{⌢}{0}}$ the n-1 zero vector and $\vec{\vec{\overset{⌢}{a}}}$ a general, symmetric $(n - 1) \times (n - 1)$ matrix. It can easily be confirmed that $\vec{\vec{a}}$ constructed as indicated by Eq. (B5) duly satisfies the conditions of Eq. (B4).

Although Eq. (B5) is helpful in identifying one way to construct $\vec{\vec{a}}$ as well as in constructively showing that the actual space where $\vec{\vec{a}}$ belongs has n(n-1)/2 dimensions [the number of independent components of the general, symmetric $(n - 1) \times (n - 1)$ matrix $\vec{\vec{\overset{⌢}{a}}}$ ] 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 $[{\vec{q}}^{(1)}, {\vec{q}}^{(2)}, \dots, {\vec{q}}^{(n - 1)}]$ ⁠. However, once we have identified, thanks to Eq. (B5) the dimensionality of the space of $\vec{\vec{a}}$ 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 $α_{i j}$ coefficients $(i = 1, 2, \dots, n, j = 1, 2, \dots, n, j \neq i)$ n × n elementary matrices ${\vec{\vec{a}}}_{i j}$ defined through the specification of their k, l component, ${({\vec{\vec{a}}}_{i j})}_{k l}$ ⁠, $(k = 1, 2, \dots, n, l = 1, 2, \dots, n)$ as

{({\vec{\vec{a}}}_{i j})}_{k l} = (δ_{k i} - δ_{k j}) (δ_{l i} - δ_{l j}) .

(B6)

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) \neq (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 $\vec{\vec{a}}$ can be described, as

\vec{\vec{a}} = \sum_{i = 1}^{n} \sum_{\begin{matrix} j = 1 \\ j \neq i \end{matrix}}^{n} α_{i j} {\vec{\vec{a}}}_{i j},

(B7)

with $α_{i j} = α_{j i}$ ⁠. If Eq. (B7) is then used to substitute for the nxn matrix $\vec{\vec{a}}$ in the general expression for the primitive symmetric dissipation bracket, Eq. (B1) we do indeed get the proposed expression, Eq. (15) which is therefore as general as the former one. QED.

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 $\int_{- t}^{x} I_{0} (\hat{z}) d \hat{x}$ ⁠, the expansion can be expressed as follows:

\begin{matrix} \int_{- t^{*}}^{x^{*}} I_{0} (\hat{z}) d \hat{x} = \int_{- t^{*}}^{x^{*}} \sum_{k = 0}^{\infty} \frac{{\hat{z}}^{2 k}}{2^{2 k} {(k!)}^{2}} d \hat{x} \\ = \sum_{m = 0}^{\infty} (\frac{(x^{*}^{2 m + 1} + t^{*}^{2 m + 1})}{16^{m} (2 m + 1) m!} (\frac{{(\frac{t^{*}}{2})}^{- m}}{2^{- m}} \sum_{n = 0}^{\infty} \frac{{(\frac{t^{*}}{2})}^{2 n + m}}{2^{2 n + m} n! (n + m)!})) \\ = \sum_{m = 0}^{\infty} (\frac{(x^{*}^{2 m + 1} + t^{*}^{2 m + 1})}{4^{m} (2 m + 1) m!} (\frac{I_{m} (\frac{t^{*}}{2})}{t^{*}^{m}})) . \end{matrix}

(C1)

Similarly, the second integral term can be expressed as

\begin{matrix} \int_{- t^{*}}^{x^{*}} \frac{I_{1} (\hat{z})}{\hat{z}} d \hat{x} = \int_{- t^{*}}^{x^{*}} \sum_{k = 0}^{\infty} \frac{{\hat{z}}^{2 k}}{2^{2 k + 1} (k!) ((k + 1)!)} d \hat{x} \\ = \frac{1}{2} \sum_{m = 0}^{\infty} (\frac{(x^{*}^{2 m + 1} + t^{*}^{2 m + 1})}{16^{m} (2 m + 1) m!} (\frac{{(\frac{t^{*}}{2})}^{- (m + 1)}}{2^{- (m + 1)}} \sum_{n = 0}^{\infty} \frac{{(\frac{t^{*}}{2})}^{2 n + m + 1}}{2^{2 n + m + 1} n! (n + m + 1)!})) \\ = 2 \sum_{m = 0}^{\infty} (\frac{(x^{*}^{2 m + 1} + t^{*}^{2 m + 1})}{4^{m} (2 m + 1) m!} (\frac{I_{m + 1} (\frac{t^{*}}{2})}{t^{*}^{m + 1}})) . \end{matrix}

(C2)

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

\begin{matrix} w_{A} (x^{*}, t^{*}) = e^{- \frac{t^{*}}{2}} + \frac{1}{2} e^{- \frac{t^{*}}{2}} [\sum_{m = 0}^{\infty} (\frac{(x^{*}^{2 m + 1} + t^{*}^{2 m + 1})}{4^{m} (2 m + 1) m!} (\frac{I_{m} (\frac{t^{*}}{2})}{t^{*}^{m}}))] \\ + \frac{1}{2} t^{*} e^{- \frac{t^{*}}{2}} [\sum_{m = 0}^{\infty} (\frac{(x^{*}^{2 m + 1} + t^{*}^{2 m + 1})}{4^{m} (2 m + 1) m!} (\frac{I_{m + 1} (\frac{t^{*}}{2})}{t^{*}^{m + 1}}))] . \end{matrix}

(C3)

References

1.

A. E.

Fick

, “

Ueber diffusion

,”

Poggendorff's Ann. Phys.

170

,

59

–

86

(

1855

), Chap. IV.

https://doi.org/10.1002/andp.18551700105

Google Scholar

Crossref

2.

L.

Onsager

, “

Theories and problems of liquid diffusion

,”

Ann. N. Y. Acad. Sci.

46

,

241

–

265

(

1945

).

https://doi.org/10.1111/j.1749-6632.1945.tb36170.x

Google Scholar

Crossref

PubMed

3.

E. L.

Cussler

,

Diffusion

:

Mass Transfer in Fluid Systems

, 2nd ed. (

Cambridge University Press

,

Cambridge

,

1997

).

Google Scholar

Crossref

4.

L.

Onsager

, “

Reciprocal relations in irreversible processes. I

,”

Phys. Rev.

37

,

405

–

426

(

1931

).

https://doi.org/10.1103/PhysRev.37.405

Google Scholar

Crossref

5.

L.

Onsager

, “

Reciprocal relations in irreversible processes. II

,”

Phys. Rev.

38

,

2265

–

2279

(

1931

).

https://doi.org/10.1103/PhysRev.38.2265

Google Scholar

Crossref

6.

H. B. G.

Casimir

, “

On Onsager's principle of microscopic reversibility

,”

Rev. Mod. Phys.

17

,

343

–

350

(

1945

).

https://doi.org/10.1103/RevModPhys.17.343

Google Scholar

Crossref

7.

J. C.

Maxwell

, “

Illustrations of the dynamical theory of gases. I. On the motions and collisions of perfectly elastic spheres

,”

Philos. Mag.

19

,

19

–

32

(

1860

).

https://doi.org/10.1080/14786446008642818

Google Scholar

Crossref

8.

J. C.

Maxwell

, “

Illustrations of the dynamical theory of gases. II. On the process of diffusion of two or more kinds of moving particles among one another

,”

Philos. Mag.

20

,

21

–

37

(

1860

).

https://doi.org/10.1080/14786446008642902

Google Scholar

Crossref

9.

J.

Stefan

, “

On the equilibrium and movement of gas mixtures, in particular diffusion

,”

Sitzungsber. Kais. Akad. Wiss. Wien

63

,

63

–

124

(

1871

).

Google Scholar

10.

J.

Stefan

, “

On the dynamical theory of diffusion of gases

,”

Sitzungsber. Kais. Akad. Wiss. Wien

65

,

323

–

363

(

1872

).

Google Scholar

11.

C. F.

Curtiss

and

J. O.

Hirschfelder

, “

Transport properties of multicomponent gas mixtures

,”

J. Chem. Phys.

17

,

550

–

555

(

1949

).

https://doi.org/10.1063/1.1747319

Google Scholar

Crossref

12.

J. O.

Hirschfelder

,

C. F.

Curtiss

, and

R. B.

Bird

, “

Molecular theory of gases and liquids

,”

Corrected Printing with Notes Added

(

John Wiley and Sons

,

New York

,

1964

).

Google Scholar

Crossref

13.

C. F.

Curtiss

and

R. B.

Bird

, “

Multicomponent diffusion

,”

Ind. Eng. Chem. Res.

38

,

2515

–

2522

(

1999

);

https://doi.org/10.1021/ie9901123

Crossref

Google Scholar

C. F.

Curtiss

and

R. B.

Bird

,

Errata

40

,

1791

(

2001

).

https://doi.org/10.1021/ie011500j

Crossref

14.

R. B.

Bird

and

D. J.

Klingenberg

, “

Multicomponent diffusion—A brief review

,”

Adv. Water Resour.

62

,

238

–

242

(

2013

).

https://doi.org/10.1016/j.advwatres.2013.05.010

Google Scholar

Crossref

15.

R. E.

Swaney

and

R. B.

Bird

, “

Transport phenomena and thermodynamics: Multicomponent mixtures

,”

Phys. Fluids

31

,

021202

(

2019

).

https://doi.org/10.1063/1.5048320

Google Scholar

Crossref

16.

C.

Truesdell

, “

Mechanical basis of diffusion

,”

J. Chem. Phys.

37

,

2336

–

2344

(

1962

).

https://doi.org/10.1063/1.1733007

Google Scholar

Crossref

17.

C.

Truesdell

and

W.

Noll

,

The Non-Linear Field Theories of Mechanics

, 3rd ed. (

Springer-Verlag

,

Berlin

,

2004

).

Google Scholar

Crossref

18.

S. H.

Lam

, “

Multicomponent diffusion revisited

,”

Phys. Fluids

18

,

073101

(

2006

).

https://doi.org/10.1063/1.2221312

Google Scholar

Crossref

19.

R.

Taylor

and

R.

Krishna

,

Multicomponent Mass Transfer

(

Wiley

,

New York

,

1993

).

Google Scholar

20.

G. D. C.

Kuiken

,

Thermodynamics of Irreversible Processes: Applications to Diffusion and Rheology

(

John Wiley and Sons

,

Chichester

,

1994

).

Google Scholar

21.

R. B.

Bird

,

W. E.

Stewart

, and

E. N.

Lightfoot

,

Transport Phenomena

, revised 2nd ed. (

John Wiley and Sons

,

New York

,

2006

).

Google Scholar

22.

R.

Datta

and

S. A.

Vilekar

, “

The continuum mechanical theory of multicomponent diffusion in fluid mixtures

,”

Chem. Eng. Sci.

65

,

5976

–

5989

(

2010

).

https://doi.org/10.1016/j.ces.2010.08.022

Google Scholar

Crossref

23.

M.

Pavelka

,

F.

Maršik

, and

V.

Klika

, “

Consistent theory of mixtures on different levels of description

,”

Int. J. Eng. Sci.

78

,

192

–

217

(

2014

).

https://doi.org/10.1016/j.ijengsci.2014.02.003

Google Scholar

Crossref

24.

A.

Morro

, “

Modelling of thermal diffusion and thermodynamic restrictions

,”

Int. J. Eng. Sci.

85

,

125

–

135

(

2014

).

https://doi.org/10.1016/j.ijengsci.2014.08.004

Google Scholar

Crossref

25.

S.

Simić

,

M.

Pavić-Čolić

, and

D.

Madjarević

, “

Non-equilibrium mixtures of gases: Modelling and computation

,”

Riv. Mat. Univ. Parma

6

,

135

–

214

(

2015

).

Google Scholar

26.

A. C.

Faliagas

, “

Mixed weak-perturbative solution method for Maxwell's equations of diffusion with Müller's partial stress tensor in the low velocity limit

,”

J. Comput. Phys.

308

,

322

–

346

(

2016

).

https://doi.org/10.1016/j.jcp.2015.12.046

Google Scholar

Crossref

27.

J.

Francu

and

J.

Mikyska

, “

An alternative model of multicomponent diffusion based on a combination of the Maxwell-Stefan theory and continuum mechanics

,”

J. Comput. Phys.

400

,

108962

(

2020

).

https://doi.org/10.1016/j.jcp.2019.108962

Google Scholar

Crossref

28.

V. G.

Mavrantzas

and

A. N.

Beris

, “

Modelling of the rheology and flow-induced concentration changes in polymer solutions

,”

Phys. Rev. Lett.

69

,

273

–

276

(

1992

);

https://doi.org/10.1103/PhysRevLett.69.273

Crossref

[PubMed]

Google Scholar

V. G.

Mavrantzas

and

A. N.

Beris

,

Errata

70

,

2659

(

1993

).

https://doi.org/10.1103/PhysRevLett.70.2659

Crossref

29.

H.

Ji

and

E.

Helfand

, “

Concentration fluctuations in sheared polymer solutions

,”

Macromolecules

28

,

3869

–

3880

(

1995

).

https://doi.org/10.1021/ma00115a017

Google Scholar

Crossref

30.

M.

Cromer

,

M. C.

Villet

,

G. H.

Fredrickson

,

L. G.

Leal

,

R.

Stepanyan

, and

J. H.

Bulters

, “

Concentration fluctuations in polymer solutions under extensional flow

,”

J. Rheol.

57

,

1211

–

1235

(

2013

).

https://doi.org/10.1122/1.4808411

Google Scholar

Crossref

31.

S.

Hooshyar

and

N.

Germann

, “

A thermodynamic study of shear banding in polymer solutions

,”

Phys. Fluids

28

,

063104

(

2016

).

https://doi.org/10.1063/1.4953859

Google Scholar

Crossref

32.

F.

Municchi

,

P.

Nagrani

, and

I. C.

Christov

, “

A two-fluid model for numerical simulation of shear-dominated suspension flows

,”

Int. J. Multiphase Flow

120

,

103079

(

2019

).

https://doi.org/10.1016/j.ijmultiphaseflow.2019.07.015

Google Scholar

Crossref

33.

N.

Germann

,

L. P.

Cook

, and

A. N.

Beris

, “

Nonequilibrium thermodynamic modeling of the structure and rheology of concentrated wormlike micellar solutions

,”

J. Non-Newtonian Fluid Mech.

196

,

51

–

57

(

2013

).

https://doi.org/10.1016/j.jnnfm.2012.12.010

Google Scholar

Crossref

34.

N.

Germann

,

L. P.

Cook

, and

A. N.

Beris

, “

Investigation of the inhomogeneous shear flow of a wormlike micellar solution using a thermodynamically consistent model

,”

J. Non-Newtonian Fluid Mech.

207

,

21

–

31

(

2014

).

https://doi.org/10.1016/j.jnnfm.2014.02.005

Google Scholar

Crossref

35.

N.

Germann

,

L. P.

Cook

, and

A. N.

Beris

, “

A differential velocities based study of diffusion effects in shear banding micellar solutions

,”

J. Non-Newtonian Fluid Mech.

232

,

43

–

54

(

2016

).

https://doi.org/10.1016/j.jnnfm.2016.03.011

Google Scholar

Crossref

36.

R. F.

Snider

, “

Irreversible thermodynamics of multicomponent fluids and its statistical mechanics basis

,”

Phys. Rev. E

103

,

032121

(

2021

).

https://doi.org/10.1103/PhysRevE.103.032121

Google Scholar

Crossref

PubMed

37.

H.

Gomez

,

I.

Colominas

,

F.

Navarrina

,

J.

París

, and

M.

Casteleiro

, “

A hyperbolic theory for advection-diffusion problems: Mathematical foundations and numerical modeling

,”

Arch. Comput. Methods Eng.

17

,

191

–

211

(

2010

).

https://doi.org/10.1007/s11831-010-9042-5

Google Scholar

Crossref

38.

D. D.

Joseph

and

L.

Preziosi

, “

Heat waves

,”

Rev. Mod. Phys.

61

,

41

–

73

(

1989

).

https://doi.org/10.1103/RevModPhys.61.41

Google Scholar

Crossref

39.

D. D.

Joseph

and

L.

Preziosi

, “

Addendum to the paper ‘Heat waves

’,”

Rev. Mod. Phys.

62

,

375

–

391

(

1990

).

https://doi.org/10.1103/RevModPhys.62.375

Google Scholar

Crossref

40.

N. S.

Kalospiros

,

B. J.

Edwards

, and

A. N.

Beris

, “

Internal variables for relaxation phenomena in heat and mass transfer

,”

Int. J. Heat Mass Transfer

36

,

1191

–

1200

(

1993

).

https://doi.org/10.1016/S0017-9310(05)80089-4

Google Scholar

Crossref

41.

A. N.

Beris

and

B. J.

Edwards

,

Thermodynamics of Flowing Systems with Internal Microstructure

(

Oxford University Press

,

New York

,

1994

).

Google Scholar

Crossref

42.

H. C.

Öttinger

,

Beyond Equilibrium Thermodynamics

(

Wiley-Interscience

,

Hoboken, New Jersey

,

2005

).

Google Scholar

Crossref

43.

B. J.

Edwards

, “

An analysis of single and double generator thermodynamic formalisms for the macroscopic description of complex fluids

,”

J. Non-Equilib. Thermodyn.

23

,

301

–

333

(

1998

).

https://doi.org/10.1515/jnet.1998.23.4.301

Google Scholar

Crossref

44.

H.

Struchtrup

and

H. C.

Öttinger

, “

Thermodynamically admissible 13-moment equations

,”

Phys. Fluids

34

,

017105

(

2022

).

https://doi.org/10.1063/5.0078780

Google Scholar

Crossref

45.

M.

Pavelka

,

V.

Klika

, and

M.

Grmela

,

Multiscale Thermodynamics

(

De Gruyter

,

Berlin

,

2018

).

Google Scholar

46.

A. N.

Beris

, “

Continuum mechanics modeling of complex fluid systems following Oldroyd's seminal 1950 work

,”

J. Non-Newtonian Fluid Mech.

298

,

104677

(

2021

).

https://doi.org/10.1016/j.jnnfm.2021.104677

Google Scholar

Crossref

47.

H.

Goldstein

,

Classical Mechanics

, 2nd ed. (

Addison-Wesley

,

Reading, MA

,

1980

).

Google Scholar

48.

L. D.

Landau

and

E. M.

Lifshitz

,

Mechanics

, 3rd ed. (

Pergamon Press

,

Oxford

,

1976

).

Google Scholar

49.

A.

Clebsch

, “

Über die Integration der Hydrodynamische Gleichungen

,”

J. Reine Angew. Math.

56

,

1

–

10

(

1895

).

Google Scholar

50.

V. I.

Arnold

, “

Sur la géometrie différentielle des groupes de Lie de dimension infini et ses applications dans l'hydrodynamique des fluides parfaits

,”

Ann. Inst. Fourier

16

,

319

–

361

(

1966

).

https://doi.org/10.5802/aif.233

Google Scholar

Crossref

51.

A. N.

Kaufman

, “

Dissipative Hamiltonian systems: A unifying principle

,”

Phys. Lett. A

100

,

419

–

422

(

1984

).

https://doi.org/10.1016/0375-9601(84)90634-0

Google Scholar

Crossref

52.

P. J.

Morrison

, “

Bracket formulation for irreversible classical fields

,”

Phys. Lett. A

100

,

423

–

427

(

1984

).

https://doi.org/10.1016/0375-9601(84)90635-2

Google Scholar

Crossref

53.

M.

Grmela

, “

Bracket formulation of dissipative fluid mechanics equations

,”

Phys. Lett. A

102

,

355

–

358

(

1984

).

https://doi.org/10.1016/0375-9601(84)90297-4

Google Scholar

Crossref

54.

B. D.

Coleman

,

M.

Fabrizio

, and

D. R.

Owen

, “

On the thermodynamics of second sound in dielectric crystals

,”

Arch. Ration. Mech. Anal.

80

,

135

–

158

(

1982

).

https://doi.org/10.1007/BF00250739

Google Scholar

Crossref

55.

S.

Gupta

,

J. D.

Schieber

, and

D.

Venerus

, “

Anisotropic thermal conduction in polymer melts in uniaxial elongation flows

,”

J. Rheol.

57

,

427

–

439

(

2013

).

https://doi.org/10.1122/1.4776237

Google Scholar

Crossref

56.

K.

Murota

,

Discrete Convex Analysis

(

SIAM

,

Philadelphia

,

2003

).

Google Scholar

Crossref

57.

T.

Matolcsi

and

P.

Ván

, “

On the objectivity of time derivatives

,”

Atti Accad. Peloritana Pericolanti Cl. Sci. Fis. Mat. Nat.

86

(

Suppl. 1

),

C1S0801015

(

2008

).

https://doi.org/10.1478/C1S0801015

Google Scholar

58.

J. G.

Oldroyd

, “

On the formulation of rheological equations of state

,”

Proc. R. Soc. London, Ser. A

200

,

523

–

541

(

1950

).

https://doi.org/10.1098/rspa.1950.0035

Google Scholar

Crossref

59.

S. I.

Sandler

and

J. S.

Dahler

, “

Nonstationary diffusion

,”

Phys. Fluids

7

,

1743

–

1746

(

1964

).

https://doi.org/10.1063/1.2746771

Google Scholar

Crossref

60.

R. J.

Bearman

and

J. G.

Kirkwood

, “

Statistical mechanics of transport processes. XI. Equations of transport in multicomponent systems

,”

J. Chem. Phys.

28

,

136

–

145

(

1958

).

https://doi.org/10.1063/1.1744056

Google Scholar

Crossref

61.

D.

Kondepudi

and

I.

Prigogine

,

Modern Thermodynamics

(

Wiley

,

Chichester

,

1998

).

Google Scholar

62.

W. M.

Deen

,

Analysis of Transport Phenomena

, 2nd ed. (

Oxford University Press

,

New York

,

2012

).

Google Scholar

63.

P. M.

Mwasame

,

N. J.

Wagner

, and

A. N.

Beris

, “

On the macroscopic modeling of dilute emulsions under flow in the presence of particle inertia

,”

Phys. Fluids

30

,

030704

(

2018

).

https://doi.org/10.1063/1.5001477

Google Scholar

Crossref

64.

P. M.

Mwasame

,

N. J.

Wagner

, and

A. N.

Beris

, “

Micro-inertia effects in material flow

,”

J. Non-Equilib. Thermodyn.

44

,

235

–

246

(

2019

).

https://doi.org/10.1515/jnet-2018-0072

Google Scholar

Crossref

65.

D. C.

Venerus

and

H. C.

Öttinger

,

A Modern Course in Transport Phenomena

(

Cambridge University Press

,

Cambridge, UK

,

2018

).

Google Scholar

66.

J. C.

Maxwell

, “

On the dynamical theory of gases

,”

Philos. Trans. R. Soc. A

157

,

49

–

88

(

1867

).

Google Scholar

67.

C.

Cattaneo

, “

Sur une forme de l'équation de la chaleur éliminant le paradoxe d'une propagation instantanée

,”

C.R. Acad. Sci. Paris

247

,

431

–

433

(

1958

).

Google Scholar

68.

V.

Ciancio

and

A.

Palumbo

, “

A thermodynamic theory with hidden vectorial variables on possible interactions among heat conduction, diffusion phenomena, viscous flow and chemical reactions in fluid mixtures

,”

Atti Accad. Peloritana Pericolanti

97

(

S1

),

A4

(

2019

).

https://doi.org/10.1478/AAPP.97S1A4

Google Scholar

69.

A. K.

Das

, “

A non-Fickian diffusion equation

,”

J. Appl. Phys.

70

,

1355

–

1358

(

1991

).

https://doi.org/10.1063/1.349592

Google Scholar

Crossref

70.

S.

Sieniutycz

and

R. S.

Berry

, “

Least-entropy generation: Variational principle of Onsager's type for transient hyperbolic heat and mass transfer

,”

Phys. Rev. A

46

,

6359

–

6370

(

1992

).

https://doi.org/10.1103/PhysRevA.46.6359

Google Scholar

Crossref

PubMed

71.

S.

Sieniutycz

, “

The variational approach to nonstationary Brownian and molecular diffusion described by wave equations

,”

Chem. Eng. Sci.

39

,

71

–

80

(

1984

).

https://doi.org/10.1016/0009-2509(84)80131-1

Google Scholar

Crossref

72.

H.

Gomez

, “

A hyperbolic formulation for convective diffusive problems in CFD

,” Ph.D. dissertation (University of A Coruña, 2006).

73.

R. C.

Mittal

and

S.

Dahiya

, “

Numerical simulation on hyperbolic diffusion equations using modified cubic B-spline differential quadrature methods

,”

Comput. Math. Appl.

70

,

737

–

749

(

2015

).

https://doi.org/10.1016/j.camwa.2015.04.022

Google Scholar

Crossref

74.

R.

Sureshkumar

and

A. N.

Beris

, “

Effect of artificial stress diffusivity on the stability of numerical calculations and the flow dynamics of time-dependent viscoelastic flows

,”

J. Non-Newtonian Fluid Mech.

60

,

53

–

80

(

1995

).

https://doi.org/10.1016/0377-0257(95)01377-8

Google Scholar

Crossref

75.

K. D.

Housiadas

and

A. N.

Beris

, “

Polymer-induced drag reduction: Effects of the variations in elasticity and inertia in turbulent viscoelastic channel flow

,”

Phys. Fluids

15

,

2369

–

2384

(

2003

).

https://doi.org/10.1063/1.1589484

Google Scholar

Crossref

2022

Author(s)

Flux-based modeling of heat and mass transfer in multicomponent systems

I. INTRODUCTION

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

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

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

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

NOMENCLATURE

Description

Greek symbols

Brackets/scripts

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

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

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

References

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

Flux-based modeling of heat and mass transfer in multicomponent systems

I. INTRODUCTION

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

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

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

V. CONCLUSIONS

ACKNOWLEDGMENTS

AUTHOR DECLARATIONS

Conflict of Interest

DATA AVAILABILITY

NOMENCLATURE

Description

Greek symbols

Brackets/scripts

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

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

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

References

Citing articles via

Resources

Explore

pubs.aip.org

Connect with AIP Publishing

This Feature Is Available To Subscribers Only