Revised Formulation of Fick's, Fourier's, and Newton's Laws for Spatially Varying Linear Transport Coefficients.

We argue that for situations involving spatially varying linear transport coefficients (diffusivities, thermal conductivities, and viscosities), the original Fick's, Fourier's, and Newton's law equations should be modified to place the diffusivity, thermal conductivity, and viscosity inside the derivative operator, that is, in one-dimensional rectilinear situations, , , and . We present simple derivations of these formulas in which diffusive mass transfer, conductive heat transfer, and molecular momentum transfer processes are described using lattice random walk models. We also present simple examples demonstrating how these modifications affect calculations.

Introduction

Adolph Fick proposed in 1855[1] that in a simple one-dimensional situation, the diffusive flux, say of species A, is proportional to the concentration gradient of the species along the system axis (the x-axis) (“Fick’s first law”)The proportionality factor D is what is called the diffusion coefficient or diffusivity. From this relationship, Fick also computed the rate of change of the concentration of A by diffusion (“Fick’s second law” or “equation of continuity for A”)where D is assumed to be spatially “constant” and “dependent (only) upon the nature of the substances”. In general situations, however, the diffusivity is position-dependent because of the spatial variation of, for instance, concentration or temperature and the dependence of the diffusivity on these variables. A common approach to deal with spatially varying diffusivity is to use the following modification of eq abovewhere use is still made of Fick’s first law equation in its original form (eq ); the only modification is that D is not factored out of the outer derivative.[2] Recently, we proposed an argument that to exactly track spatial changes in D, Fick’s first law equation must also be generalized to the form[3]which thus also gives a different equation for Fick’s second lawWhen mass flows have components in all directions, eqs and 5 should be expressed in the vectorial form, respectively, asandThe correctness of eq (or eq ) is intuitively obvious; eq predicts that even in the absence of a concentration gradient, a diffusive flux of A arises if there exists a gradient of diffusivity and also that in the absence of a net material flux (j = 0), a spatial gradient of diffusivity results in a spontaneous buildup of nonuniform concentrations of A (note that these phenomena are not describable within the original Fick’s law framework).[3]

As discussed in ref (3), the above eq can be obtained using the lattice model description of diffusion.[4] For the purpose of setting the basis for the discussion to be presented in the present article, this derivation is briefly repeated here. In this derivation, we consider a steady-state situation described in Figure A, where a suspension of Brownian particles (molecules) (species A) is contained in a tube. The concentration of this species varies only along the tube axis (x-direction). For simplicity, the three-dimensional (3D) space within the tube is pictured as being divided into square lattice sites; each site of the lattice can hold at most one molecule at a time and has a characteristic dimension of λ. Under this setting, we model the diffusive motion of the molecules to be a 3D random walk process.[4] Specifically, we assume that each molecule steps to an adjacent lattice site with a jump frequency of ν. As illustrated in Figure A, let us now consider a plane of constant x between two adjacent lattice layers. The magnitude of the flux of the molecules from lesser x to greater x across the plane (located at x) and that of the flux from greater x to lesser x across the same plane can be calculated, respectively, asandwhere c is the concentration of A (e.g., in units of mass of molecules per unit volume), and the one-sixth factor is introduced because in 3D space, a molecule can move to one of the six nearest-neighbor sites with the equal probability of 1/6. Therefore, the net mass flux of the molecules across the plane located at x is given bywhere the product is the diffusion coefficient (in lattice units) for that location, i.e., . The resulting equation (eq ) is the steady-state version of eq . This derivation has been given in an earlier publication (in the Supporting Information of ref (3)).

Figure 1

Molecular transport of (A) mass, internal energy, and (B) y momentum due to gradients of Dc, kT, and μv along the x-direction, respectively.

The realization of this generalized form of Fick’s first law equation raises two important questions. First, can (or even should) this same argument be applied to Fourier’s law of heat conduction and Newton’s law of viscosity, respectively, for situations involving spatially varying thermal conductivities and viscosities? The answer is intuitively obvious (yes) because as stated in Bird, Stewart, Lightfoot, and Klingenberg (BSLK),[2] “the molecular mechanisms responsible for the transport of chemical species, energy and momentum are closely related (the same molecular motions and interactions are responsible for diffusivity, thermal conductivity and viscosity)”. The second question is: In reality, how much difference would using, for instance, eq instead of eq (for solving mass transfer problems involving position-dependent diffusivities) make to the result? Or, alternatively put, what is the range in which the commonly used approximation (eq , as opposed to the more accurate equation, eq ) is valid? In the context of these questions, the present paper attempts to serve two purposes. It first presents simplistic arguments that justify the generalization of Fourier’s thermal conductivity and Newton’s viscosity equations, respectively, for spatially varying thermal conductivities and viscosities; these arguments are similar to those used above for generalizing Fick’s first law equation for spatially varying diffusivities. It then discusses simple (but realistic) examples demonstrating how these modifications impact the calculations and predictions of the equations of change for mass, energy, and momentum.

Results and Discussion

Derivations

Fick deduced his first law of diffusion (eq ) by analogy with Fourier’s law of heat conduction (and Ohm’s law of electrical conduction).[1] Likewise, it is reasonable to expect that the same generalization as in eq is applicable to Fourier’s law. To actually show such a derivation, let us first consider heat conduction in a gaseous system. We use the same lattice model description of diffusion as in Section ; see Figure A for the geometry of the system. Assuming that the lattice dimension (λ) is comparable to the mean free path of the molecules, the magnitude of the heat flux (i.e., the molecular transport of the internal energy due to collision of the molecules) from lesser x to greater x across the plane located at x and that of the heat flux from greater x to lesser x across the same plane can be calculated, respectively, asandwhere C̅ is the heat capacity per molecule (at a constant volume) and Tref is the reference temperature; C̅(T – Tref) thus gives the internal energy per molecule. All other parameters are the same as in eqs and 9. The net heat flux across the plane located at x is, therefore, given bywhere is the thermal conductivity; Tref is set to 0. Note that the exactly same argument can also be applied to heat conduction in a solid body, simply by replacing c with the phonon concentration, νλ with the mean phonon velocity, and C̅ with the phonon heat capacity.[5] Therefore, when the thermal conductivity varies with position, the original Fourier’s law equationshould be generalized to the formwhich yields a different expression, for instance, for the equation of temperature for a solidsee eqs (11.2–10) of ref (2) for the original version of the equation. When heat flows have components in all directions, eq should be expressed in the vectorial form aswhere a bold character is used to denote a vector.

It is trivial to show the same derivation for momentum transfer in a lattice gas (flowing in the y-direction with a velocity gradient dv/dx). Again assuming that the lattice dimension (λ) is comparable to the mean free path of the molecules, the magnitude of the y momentum flux (i.e., the molecular transport of the y momentum due to collision of the molecules) across the plane of constant x located at x in the positive x-direction and that of the y momentum flux across the same plane in the negative x-direction (Figure B) can be calculated, respectively, asandwhere m is the mass of the molecule and v is the velocity of the gas along the y-direction; mv thus gives the y momentum of the molecule. The net y momentum flux across the plane of constant x located at x is, therefore, given bywhere is the (shear) viscosity of the gas. Sir Isaac Newton proposed in “Principia” that “the resistance (τ) arising from the lack of slipperiness (μ) in a fluid is proportional to the velocity with which the parts of the fluid are separated from one another (∂v/∂x)”,[6] which has been formulated later by scientists into the equationFor spatially varying viscosities, this original Newton’s viscosity equation is generalized to the formLogic suggests that for flows involving position-dependent viscosities the full vector–tensor expression for the viscous stress (momentum flux) tensor should also be generalized to the formwhere [∇(μv)]† is the transpose of the ∇(μv) tensor, κ is the dilatational viscosity, δ is the unit tensor, and bold characters denote vector and tensor quantities; further study is needed to prove this generalization rigorously. Accordingly, the equation of motion[2] becomeswhere ρ is the density, p is the pressure, g is the gravitational acceleration, and D/Dt is the substantial time derivative operator. Note that now even for an incompressible fluid (∇·v = 0) the third and fourth terms on the right-hand side above do not vanish in general, which will make computation more difficult.

Implications

To our knowledge, the generalized expressions of Fick’s law of diffusion, Fourier’s law of heat conduction, and Newton’s law of viscosity proposed in the present work (eqs , 17, and 23, respectively) have not been demonstrated in the transport phenomena/continuum mechanics literature previously (although such formalisms have been implied in statistical mechanics as briefly discussed in the next subsection). For instance, in COMSOL (a commercial finite element method simulator that is widely used for solving fluid mechanics problems), nonisothermal flow problems (involving spatially varying μ and k) are typically solved using the modified equations of motion and temperature, in which (analogously to the diffusive flux term in eq , that is, −∇·j = ∇·(D∇c) in full vector notation) the viscous momentum flux term (−∇·τ) and the conductive heat flux term (−∇·q) are calculated, respectively, using the original forms of Newton’s law of viscosity and Fourier’s law of heat conduction[7]where the simplicity of the latter expression is due to the assumption of constant ρ (∇·v = 0), andWe now argue that these conventional expressions (eqs and 26) are only approximations for the general expressions given in eqs and 17, respectively. It will require extensive investigations to establish the ranges of conditions under which the use of the generalized formulas that we propose (eqs and 17) is required rather than the standard “approximations” (eqs and 26). In the present paper, we do not intend to offer a full analysis of this question. Instead, we will present simple examples that demonstrate a need for further research in this direction. Ordinary examples are well suited for this purpose. For this reason, examples have been taken from one of the most popular textbooks of transport phenomena, BSLK.

The first example concerns the heating of electric wire (Figure A) with temperature-dependent thermal and electrical conductivities, k and ke, respectively (Problem #10C.1 of BSLK)where k0 and ke0 are, respectively, the thermal and electrical conductivities at a reference temperature T0, Θ(=(T – T0)/T0) is a dimensionless temperature, and the coefficients α and β (i = 1, 2,...) are material-dependent constants. The steady-state energy balance in cylindrical coordinates gives[2]where E is the voltage drop and L is the wire length; note that the quantity on the right-hand side of the equation represents the rate of heat generation per unit length of the wire due to electrical energy dissipation, and the final expression is obtained by substituting the original Fourier’s law equation (similarly to eq , in cylindrical coordinates) for q. When this equation is solved using a perturbation method with the boundary conditions that T is finite at r = 0 and T = T0 at r = R, one obtains a solution for the radial temperature profile in the electrically heated wire in a dimensionless formwhere B = (ke0R2E2)/(k0L2T0), ξ = r/R, and O(B2) means terms of the order of B2 and higher.[2] If we use the generalized form of Fourier’s law (eq , that is, in cylindrical coordinates), the energy balance equation (eq ) changes toNow, it can be shown that this modification leads to a (slightly) different solutionBy way of example, for a copper wire of radius R = 2 mm and length L = 5 m across a voltage drop of E = 40 V at T0 = 20 °C (k0 = 385 W/(m K), α1 = 0.0508 (estimated from Tables 9.5–4 of BSLK), ke0 = 5.99 × 105 Ω–1 cm–1, β1 = 0.872 (estimated from Tables 9.5–4 and eqs 9.9–1 of BSLK), which gives a value of B = 0.136 for the dimensionless heat source strength),[2] the temperature distributions were calculated using the two equations above, eqs and 32. These results are displayed in Figure B. As shown in the figure, for these mild parameter values used, the generalized Fourier’s law produces predictions for electrical heating of the wire that are practically indistinguishable from those of the original Fourier’s law; the differences are less than 0.2%, although this small difference increases as the rate of electrical energy dissipation (B) is increased. This result is due to the fact that, as shown in Figure C, the heat flux due to the thermal conductivity gradient is negligible relative to the heat flux due to the temperature gradient .

Figure 2

(A) Electrically heated wire. (B) Dimensionless temperature profiles predicted using eqs vs 32. Difference (%) is defined as (Θ(eq ) – Θ(eq ))/Θ(eq ) × 100. (C) Dimensionless deviation factor as a function of r/R; absolute magnitudes are compared because the and terms have opposite signs.

Next, let us consider a nonisothermal momentum transfer process that involves a (Newtonian) liquid flowing downward (in the positive y-direction) along the surface of a vertical plane in a steady laminar flow (Figure A); this example is again taken from BSLK (examples 11.4–3 and 2.2–2). The temperature of the free liquid surface (x = 0) is kept at a constant T0 and that of the solid surface (x = δ) is kept at Tδ. At these temperatures, the liquid has viscosities of μ0 and μδ, respectively. For simplicity, we assume that within this given range of temperature, the density (ρ) and thermal conductivity (k) of the liquid are constant. Due to the Arrhenius-type dependence of viscosity on temperature (that is, μ/μ0 = exp[B(1/T – 1/T0)], where B is a constant), the spatial variation of viscosity also has an exponential character[2]Substitution of the original Newton’s viscosity law (eq ) with variable viscosity (eq ) into the steady-state y momentum balance giveswhich upon integration with the no-slip boundary condition (that is, v = 0 at x = δ) gives[2]If we use the generalized Newton’s law of viscosity (eq ), the above momentum balance equation (eq ) is changed toThis equation yields a different velocity profile under the same no-slip boundary conditionNote that in the constant viscosity limit (that is, when α = 0), both eqs and 37 reduce to an identical formwhich supports consistency between the two equations. For demonstration of the difference between the predictions based on the original vs generalized Newton’s law equations, we now assume that the liquid is an oil whose viscosity is μ0 = 0.16 Pa s and density is ρ = 0.8 × 103 kg/m3 at temperature T0 = 20 °C and that the falling film has a thickness of δ = 2.5 mm, the vertical wall is kept at a temperature of Tδ = 10 °C, and the Arrhenius activation energy for viscosity has a value of B = 1.04 × 103 K (for n-heptane[8]) (which gives a value of −0.125 for the dimensionless constant α). For these parameter values, the velocity profiles were calculated using the two different versions of the velocity equation shown above. As shown in Figure B, eqs and 37 predict significantly different velocity profiles. This difference further increases as the temperature gradient (α) is increased. As shown in Figure C, near the free liquid surface (at small x), the y momentum flux due to the viscosity gradient is, in fact, comparable in magnitude to the y momentum flux due to the velocity gradient .

Figure 3

(A) Nonisothermal falling liquid film. (B) Velocity profiles predicted using eqs vs 37. Difference (%) is defined as (v(eq ) – v(eq ))/v(eq ) × 100. (C) Dimensionless deviation factor as a function of x/δ; absolute magnitudes are compared because the and terms have opposite signs.

Finally, let us discuss a mass transfer example (discussed in Section 18.3 of BSLK). As shown in Figure A, a solid sphere of potassium permanganate (KMnO4) is placed in a stationary reservoir of water. KMnO4 is only slightly soluble in water; the solubility of KMnO4 in water is about 0.0758 g cm3 at 25 °C.[9] Therefore, the outward flux of the dissolved MnO4– ions (away from the sphere surface) is predominantly diffusive (the convective flux is negligible). Also, the slowness of the dissolution process allows us to make a quasi-steady-state approximation; the mass transfer process can be approximated as occurring at a steady state. The steady-state mass balance for dissolved MnO4–1 ions in spherical coordinates can be written aswhere the total mass flux of MnO4–1 (n (=cv + j)) is equated to the diffusive mass flux of MnO4– (j) in the absence of convection (v = 0), and then, j is replaced by , as related by the original Fick’s first law equation in spherical coordinates (similarly to eq ). The diffusivity for binary liquid mixtures is typically non-negligibly dependent on species concentration. By substituting an approximate expression for the concentration dependence of the diffusivity[10]in which D0 is the diffusion coefficient in infinite dilution (=1.632 × 10–5 cm2/s for MnO4– in water at 25 °C[9]) and v̲ is the specific volume of the MnO4–1 ion (estimated to be about 0.3699 cm3/g) into eq , we obtainIntegration of this equation with the boundary conditions that c = cR (which we assume to be equal to the solubility of KMnO4 in water, that is, 0.0758 g cm3 at 25 °C) at r = R (assumed to be 1 μm) and c = 0 at r → ∞ gives the concentration profile of MnO4– at the (quasi)-steady stateAlternatively, we can use the generalized Fick’s first law (eq ) to obtainUpon integration of this equation with the same boundary conditions, we getThe two expressions for (eqs and 44) are plotted as a function of in Figure B. The two curves are supposed to have common asymptotes for small and large (as required by the boundary conditions), but in fact, they are found to be practically identical to each other at all distances. This result justifies the use of the original Fick’s law to model the mass transfer associated with dissolution and diffusion of slightly soluble material in a solvent. In this case, the mass flux due to the diffusivity gradient is negligible in magnitude compared to the mass flux due to the concentration gradient (Figure C).

Figure 4

(A) Diffusion from a slightly soluble sphere. (B) Concentration profiles predicted using eqs vs 44. Difference (%) is defined as (c(eq ) – c(eq ))/c(eq ) × 100. Note that all points and curves overlap with one another. (C) Dimensionless deviation factor as a function of r/R; absolute magnitudes are compared because the and terms have opposite signs.

Consistency

We would like to add a few remarks regarding whether the proposed correction to Fick’s, Fourier’s, and Newton’s laws fits within (or violates any of) the established principles of thermodynamics and statistical physics. First, we note that eq (or eq ) suggests an interesting possibility that the net conductive heat flow may occur even in the direction of increasing temperature if the magnitude of the T∇k term is greater than the magnitude of the k∇T term (e.g., at very high temperature); this is possible, for instance, in solids because (as can be seen from eq ) the thermal conductivity of a solid typically decreases with increasing temperature, and thus ∇T and ∇k have opposite signs. We note that this prediction, though somewhat counterintuitive, does not violate the second law of thermodynamics (“no process is possible, which consists solely in the transfer of heat from one temperature level to a higher one”[11]). The second law of thermodynamics concerns processes that start and end with equilibrium states within globally isolated systems (e.g., heat exchange between two heat reservoirs), whereas Fourier’s law is an energy balance equation for a local differential control volume, which is, by definition, a non-isolated (e.g., open) system. Therefore, it is generally impertinent to discuss whether predictions of Fourier’s law are consistent with the second law of thermodynamics.

As noted in Section (also discussed in ref (3)), eq (or eq ) implies that even in the absence of a concentration gradient (∇c = 0), the net material flow occurs when the diffusivity gradient is nonzero (j = −c∇D). This situation is not unphysical. It is the chemical potential gradient that actually drives diffusion (not the concentration gradient), and a uniform concentration does not necessarily mean that the chemical potential (μ) is uniform. Note thatwhere μo and γ are the standard state chemical potential and activity coefficient of the solute, respectively, and R is the universal gas constant. Therefore, even under constant μo, T and c, a nonzero ∇μ may develop if γ varies spatially (due to spatially varying D); it is known that γ and D are related by the following equation[12]where kB is Boltzmann’s constant and f is the friction factor of the the solute molecule. Note that this argument is relevant to non-dilute situations. However, even in the dilute limit, a non-zero ∇μ may develop if μo varies spatially, for instance, because of spatially varying solvation states of a solute molecule as in the pH-phoretic situation discussed in ref (3).

Finally, we would like to point out that the proposed correction to Fick’s first law is implied in the form of the Fokker–Planck equation derived using Ito’s stochastic calculus (also known as the Kolmogorov forward equation) for the probability density of a stochastic process;[13] in one-dimensional situationswhere p is the probability density, v is the average (convection) velocity, and D is the diffusion coefficient. Similarly to eq , a probability continuity equation can be written in terms of the probability flux (jp)Comparing eq with eq givesIn the absence of convection (v = 0), eq reduces to a form analogous to eq (generalized Fick’s first law)Interestingly, if we assume that the system obeys the principle of detailed balance (microscopic reversibility (jp = 0 at all x), which is a sufficient condition for equilibrium), we obtain[13]where pe is the probability density at equilibrium. Assuming that pe is a weak function of x, we getSubstitution of this equation into eq giveswhich coincides with the conventional form of Fick’s second law equation for position-dependent diffusivity (eq ); thus, the original Fick’s second law can be considered a restricted form of the Fokker–Planck equation. We note that detailed balance is a sufficient condition for entropy maximization in an isolated system.[14] Therefore, detailed balance is not generally necessary for the continuity (Fokker–Planck) equation to be consistent with the second law of thermodynamics.

Conclusions

The diffusion(-like) equations representing Fick’s, Fourier’s, and Newton’s laws were originally derived for constant diffusivity, thermal conductivity, and viscosity, respectively. However, even when dealing with problems involving spatially varying diffusivity, thermal conductivity, and viscosity, the original Fick’s, Fourier’s, and Newton’s laws have always been used (i.e., in their original forms) without questioning the validity of such uses. We here argue that for position-dependent diffusivity, thermal conductivity, and viscosity, Fick’s, Fourier’s, and Newton’s law formulas should, in principle, be changed such that the viscosity, thermal conductivity, and viscosity are moved inside the derivative (gradient) operator; that is, in one-dimensional situations, for instance, , , and , respectively. Our examples demonstrate that even with moderate spatial variations of diffusivity, thermal conductivity, or viscosity, the proposed modifications of Fick’s, Fourier’s, and Newton’s law equations might lead to predictions that are discernibly different from those of the original formulas under certain circumstances. This issue is expected to become more important, for instance, for highly nonisothermal processes, particularly, those that involve viscous dissipation of energy that induces large spatial temperature gradients such as in fluids around rapidly moving objects. In a previous publication, we have shown that there exists a special situation in which a spontaneous build-up of non-uniform concentrations of a solute occurs in the absence of net material flux because of a spatial gradient of diffusivity, that is, ;[3] such a phenomenon (“pH phoresis” discussed in ref (3)) cannot be described using the original Fick’s first law equation. Further study on this general topic is desirable.