Dissipation Mechanisms for Fluids and Objects in Relative Motion Described by the Navier-Stokes Equation.

This work demonstrates that an additional resistance term should be included in the Navier-Stokes equation when fluids and objects are in relative motion. This is based on an observation that the effect of the microscopic molecular random velocity component parallel to the macroscopic flow direction is neglected. The two components of the random velocity perpendicular to the local mean flow direction are accounted for by the viscous resistance, e.g., by Stokes' law for spherical objects. The relationship between the mean- and the random velocity in the longitudinal direction induces differences in molecular collision velocities and collision frequency rates on the up- and downstream surface areas of the object. This asymmetry therefore induces flow resistance and energy dissipation. The flow resistance resulting from the longitudinal momentum transfer mode is referred to as thermal resistance and is quantified by calculating the net difference in pressure up- and downstream the surface areas of a sphere using a particle velocity distribution that obeys Boltzmann's transport equation. It depends on the relative velocity between the fluid and the object, the number density and the molecular fluctuation statistics of the fluid, and the area of the object and the square root of the absolute temperature. Results show that thermal resistance is dominant compared to viscous resistance considering water and air in slow relative motion to spherical objects larger than nanometer-size at ambient temperature and pressure conditions. Including the thermal resistance term in the conventional expression for the terminal velocity of spherical objects falling through liquids, the Stokes-Einstein relationship and Darcy's law, corroborates its presence, as modified versions of these equations fit observed data much more closely than the conventional expressions. The thermal resistance term can alternatively resolve d'Alembert's paradox as a finite flow resistance is predicted at both low and high relative fluid-object velocities in the limit of vanishing fluid viscosity.

Introduction and Aim

The Navier–Stokes (N–S) equation without body forces can, for incompressible fluids where ∇·u = 0, be written as[1−5]where ρ is the fluid density, t is the time, u is the local mean fluid velocity, p is the pressure, μ is the fluid viscosity, and σ is the stress tensor. Isothermal conditions are assumed; therefore the viscosity is constant. The stress tensor elements describe flow resistance due to changes in pressure and shear stresses for a fluid in motion due to internal friction between adjacent fluid elements referred to as the transversal momentum transfer mode herein, leading to the well-known viscous resistance force. The viscous resistance force on a fluid in slow relative motion to a spherical object can be calculated analytically in the continuum limit from eq using, e.g., Stokes’ stream function[6,7] method (Appendix A). It is given aswhere ζS = 6πμR is the Stokes friction factor, R is the radius of the sphere, and u1 is the magnitude of the nonperturbed relative velocity between the fluid and the sphere. The pressure contribution to the resistance in eq is the focus herein. It results from a symmetric pressure difference between the up- and the downstream part of the sphere (Appendix A), where the meaning of up- and downstream surface areas is illustrated in Figure . The symmetric imbalance in the pressure creates a net viscous resistance caused by the macroscopic velocity field equal to 2πμRu1 along spatial direction 1, which is referred to as the symmetric pressure contribution to flow resistance. The contribution from the shear stresses amounts to 4πμRu1, giving the total viscous resistance caused by a sphere as in eq . It is important to note that the viscous resistance in eq is completely determined by the transversal momentum transfer mode associated with the macroscopic flow velocity field determined from the N–S equation.

Figure 1

Up- and downstream surface areas of a sphere in blue and red, respectively, for a fluid passing from left to right along spatial direction 1. The terminology is also used for a case when a sphere moves from right to left relative to a fluid.

There are, however, both experimental and physical arguments indicating that the flow resistance for fluids and objects in relative motion should also depend on terms other than only the viscous resistance given by Stokes’ law. It has, e.g., been shown experimentally that it is difficult to predict correct terminal velocities for varying particle sizes[8−11] when using Stokes’ law as the only resistance term. Furthermore, it is challenging to fit measured molecular diffusion coefficients when varying the system temperature using the Stokes–Einstein relationship.[12−18] The absolute permeability of porous media can be interpreted to depend on absolute temperature for small laboratory-size rock samples[19−23] and related to the utilization of geothermal energy in large-scale hot water injection operations.[24−26] The resistance terms used in all of these laws/relationships are derived from the stress tensor in the N–S equation in the limit of low relative velocities proportional to the force in eq . The two latter examples indicate that an additional term should depend explicitly on the absolute temperature. The physical argument mentioned above is related to the fluctuation-dissipation theorem[27−30] (FDT) and is discussed later.

The aim of this paper is to show that there is an additional resistance term to the viscous resistance, which should be included so that the total flow resistance when gases (limitation related to gases is later relaxed to also include liquids) and objects are in slow relative motion becomes a sum of two terms. The total velocity of individual gas molecules on the microscopic level is given by v = u + w, where u is the local mean velocity defined in eq and w is the random velocity, also called peculiar velocity,[2,3,5] associated with the thermal kinetic energy of the molecules. The additional resistance term originates from the observation that the effect of the two random velocity components perpendicular to the local mean velocity is accounted for by the transversal momentum transfer mode whereas the component parallel to the macroscopic flow direction is neglected. The latter mode is referred to as the longitudinal momentum transfer mode. The pressure and the flow resistance associated with the mode are called the asymmetric pressure contribution to flow resistance and the thermal resistance,[19] respectively. Hence, thermal resistance adds to the viscous resistance; therefore, the total flow resistance for fluids and objects in relative motion becomes a sum of two contributions: (1) viscous, caused by shear stresses plus the symmetric pressure contribution and (2) thermal, caused by the asymmetric pressure contribution.

The physical origin of the thermal resistance and its quantification are given in the next section before testing and comparing its consequences to experimental data. The laws/relationships compared comprise the conventional and modified expressions for terminal velocity of spherical particles falling through water,[8−11] the Stokes–Einstein relationship,[12−18] and Darcy’s law.[31] All of these laws/relationships apply a resistance force term derived either directly or indirectly from the N–S equation. A discussion related to d’Alembert’s paradox[32] is given before summarizing the work and drawing conclusions.

Thermal Resistance

Origin of the Thermal Resistance

The physical origin of the longitudinal mode is that the momentum transfer from individual gas molecules to an object becomes a function of position in the macroscopic flow direction, e.g., along direction 1. The presence of a mean gas velocity induces differences in molecular collision velocity and collision frequency rates on the up- and downstream surface areas[28,33,34] in the direction parallel to the macroscopic velocity direction as illustrated in Figure . The reason for this is that the collisions are more violent and frequent on the upstream surface as compared to the downstream surface.[28,33,34] Since u1 and w1 are collinear and therefore correlated, the net resisting force emerging from the longitudinal mode is expected to be dominant as compared to the transversal mode, which emerges from uncorrelated velocities. That is, u and w in the latter case are always perpendicular to each other described by the off-diagonal terms in the deviatoric stress tensor where i ≠ j. Expressions for quantifying the asymmetric pressure contribution to flow resistance are developed in the next section by applying a particle velocity distribution function, which obeys Boltzmann’s transport equation (BTE)[3,35] and therefore includes both the local mean- and the random velocities, respectively.

Figure 2

Net resisting force in macroscopic flow direction 1 on a stagnant sphere caused by differences in collision velocity and collision frequency rates on the up- and downstream surfaces of a sphere.

Microscopic Interpretation of the Pressure Terms

To interpret the terms in the stress tensor microscopically and thereby establish an expression to quantify the effect of the longitudinal momentum transfer mode, a particle velocity distribution function, fNE(v1,x1), along direction 1 obeying the BTE is applied.[2,3,35] The distribution is characterized by the local mean velocity u, which is the time average of the total velocity, u = ⟨v⟩, and the random velocity, w, which has a time average equal to zero.[2,3,5] The magnitude of the random velocity is in general assumed to be much larger than the mean velocity; therefore, |w| ≫ |u|. A typical value for |w| is the root mean square of the kinetic energy known as the thermal velocity,[3], where kB is Boltzmann’s constant, T is the absolute temperature, and m is the molecular mass. The distribution function, fNE(v1,x1), is given as[2,3,5]where w1′ is the random velocity along direction 1. The distribution is only in local equilibrium in the position space and nonisotropic in the three-dimensional velocity space. It obeys mass conservation locally after eliminating the local mean velocity, u1, by a Galilean transformation.[5] It should be mentioned that the corresponding particle velocity distribution for a gas at rest is given by .[2,3,5] The latter is in global equilibrium in position space and distributed Maxwellian in velocity space. The microscopic expression for the pressure terms in the stress tensor for a gas moving with constant velocity u1 then becomesThe pressure on a unit surface of an object due to molecular collisions perpendicular to direction 1 is then given in one dimension (1D) asBefore quantifying the thermal resistance, it is beneficial to reformulate eq to the physics in the molecular collision processes by a shift in variables. Let w1 = w1′ + u1; therefore, eq becomesIf the total velocity of a gas molecule instead had been v1 = u1 – w1′, then the resulting pressure corresponding to eq would have beenIt is shown in the next section that the expressions, left and right in eqs and 7, correspond exactly to the forces obtained by considering energy and momentum conservation in gas molecule–object collisions on the left and right side of an object. The difference between these two pressures therefore gives the net resistance force due to the longitudinal momentum transfer mode.

Quantification

The term referred to as thermal resistance herein has been quantified by several authors in various studies.[27,28,33,34] The description reported by Molina[33] is followed here, whereas the terminology applied by Grassia[28] is used to classify collision processes. Consider a simple system in a one-dimensional Cartesian coordinate system comprising a large particle with mass M surrounded by a large collection of smaller particles, all with mass m. The large particle is in relative motion along spatial direction 1, u1, to the mean velocity of the smaller particles, as shown in Figure . The area of the large particle is A⊥, which is perpendicular to the direction of the mean velocity of the small particles. The large particle can, e.g., represent a grain in a porous medium surrounded by air or a Brownian particle moving in the air. The smaller particles then represent the air’s gas molecules. The relative instantaneous velocity between the large particle and an individual gas molecule is u1 + w1. There are basically two types of small–large particle collisions that can take place called catch-up (downstream part) and head-on[28] (upstream part) collisions. A catch-up collision is illustrated in (a) where a small particle catches up with a large one due to the relative velocity (w1 – u1). A head-on collision is illustrated in (b) where the large and the small particle approach each other with relative velocity (u1 – w1). Figure represents possible collision scenarios when a large particle moves from left to right through a “sea” of smaller gas molecules or when a gas flows from right to left around a stagnant large particle. The relative velocities are the same in both cases. There are two mechanisms that contribute to a net resistance force in these two collision types.

Figure 3

Catch-up (a) and head-on (b) collision mechanisms. Positive direction to the right.

The decrease in velocity of the large particle during a head-on collision exceeds the corresponding velocity increase during a gentler catch-up collision.

The head-on collision rate frequency exceeds the corresponding catch-up collision frequency.

Following Molina,[33] considering the situation in Figure where the gas molecules are described by a Maxwellian velocity distribution function, . The average force acting on the left side of the large particle due to catch-up collisions with relative velocity, w1 – u1, in (a) iswhich is equal to eq , except for the (1 + ε) term. ε is the coefficient of restitution in the collision and included to account for potential variation in momentum transfer during collisions. ε is equal to unity for fully elastic collisions, which is always assumed herein. Similarly, the average force on the right side of the large particle due to head-on collisions with relative velocity, u1 – w1, in (b), is equal to eq augmented with the (1 + ε) term. HenceThe net resisting force, FT, acting on the large particle is the difference between FLeft and FRightwhere erf is the error function. The resisting force is negative since it points to the left. Equation is the main result herein. It quantifies the incremental flow resistance term caused by the asymmetric pressure contribution to flow resistance and is always present when u1 > 0. The square root dependency related to temperature is not dependent on the assumption of a Maxwellian velocity distribution[28] but follows from the relationship between kinetic energy and absolute temperature.[36] Since the relative velocity between the large particle and the gas is assumed to be much smaller than a typical gas molecule velocity, therefore u1 ≪ VT and eq can be approximated as[33]Equation quantifies the thermal resistance force, FT u, on the gas in the low relative motion limit. It depends on the molecular number density assumed to be constant in space, the relative velocity, the area and shape of the object, and the absolute temperature. Furthermore, thermal resistance is sensitive to the system pressure via the number density term.

The longitudinal momentum transfer mode resulting in the flow resistance given by eq has also been shown to be in line with[28] the FDT,[29,30] considering conservation of momentum and energy in collisions on the microscopic level. This is not surprising since gases obviously are among the systems demonstrating fluctuations. Such fluctuations impact larger objects (e.g., Brownian particles[30]) residing in the gas with random push forces, which can be characterized statistically.[29,30] But the same fluctuations are also responsible for the systematic resistance force the object experiences when moving relatively to the gas[3,30] given by eq . A major quantity characterizing the FDT is the absolute temperature of the system.[29,30] It can therefore be stated that the thermal resistance force corresponds to a systematic resistance force to relative motion in line with the FDT. Hence, it is always present and contributes to flow resistance whenever fluids and objects are in relative motion as thermal fluctuations are always present in fluids.

The total friction factor for a sphere is then the sum of the viscous and the thermal contributions where the latter is obtained from eq (Appendix B). HenceEquation is an analytical expression for the total friction factor acting when gases are in slow relative motion to spherical objects. Its empirical consequences can therefore be compared directly to the results obtained when only using the conventional Stokes friction factor, ζS. A challenge appears, however, when predictions from eq are to be compared to experimental data, since it is strictly only valid for slow relative velocities and low-density gases. If objects fall through gases, the low resistance induces Reynolds numbers above unity, which typically is the upper limit for Stokes flow[6] to be valid. To test the consequences of eq under valid flow conditions, it is hypothesized that it also can be used for objects in relative motion to liquids. Using eq for liquids introduces, however, more uncertainty as higher-order terms in the microscopic expression for the stress tensor terms cannot be neglected anymore because of more long-range molecular interaction forces in liquids as compared to gases.[3] It is, however, common practice to apply the N–S equation interchangeably for gases and liquids[4] (referred to as the “applied approach” herein), and the practical results obtained are quite satisfactory,[4] even though not all of the theoretical details are completely worked out on the molecular level. Hence, the hypothesis above related to the use of eq for liquids is investigated according to the current “applied approach”. That is, it is tested and compared to observed data and the conclusion about its status is based on the results, although not all details about the exact molecular level processes are worked out.[1,3] Furthermore, since the energy dissipated cannot be calculated accurately from the N–S equation ab initio due to the use of constituent relationships and approximations,[1,3] neither for the conventional nor for the new extended approach, a fit parameter is used for both cases to equate the calculated initial values to a reference measured data point. The fit parameter is higher for the extended case as the energy dissipated now is distributed into two rather than one term. Laws/relationships to be tested include the following:

terminal velocity of spherical particles[6,7] falling through liquids,

the Stokes–Einstein relationship[37] describing the diffusion of ions/molecules through solvents, and

Darcy’s law[1,31] describing low-rate single-phase fluid flow through porous media.

Additionally, d’Alembert’s paradox[32] stating that no drag force exists on an object residing in a flowing fluid with vanishing viscosity is discussed.

The assumptions related to the application of eq if not otherwise stated are (1) low relative motion between object and fluid, Reynolds numbers typically <1; (2) steady relative motion, therefore ; (3) fluid and object in thermal equilibrium and the object is neither attractive nor repelling with respect to the fluid molecules; (4) no accumulation of fluid molecules on the object surface and all molecular collisions totally elastic; (5) incompressible fluids and objects; (6) changes in fluid density and mean velocity on the macroscopic level occur slowly compared to processes on intermediate time and length scales,[2,5] hence the concept of local equilibrium is appropriate; and (7) change in internal energy of the fluid is caused by a change in translation energy only, excluding changes in other internal modes such as rotation and vibration.[3]

Consequence Testing Results

Modification of the Terminal Velocity Expression

A balance of the gravity force with the sum of buoyancy and frictional forces (eq ) according to Newton’s second law for a sphere with density ρS falling through a liquid with density ρ gives the terminal velocity uT. It is conventionally expressed aswhere βT is a fit parameter, which is related to the argument in the previous section. It also includes effects related to nonsphericity[6,8] and surface properties of the sphere, e.g., adsorbed material, irregularities, and polarity compared to the fluid molecules.[38,39]Equation shows that the terminal velocity expression in eq should be modified. The expression for the modified terminal velocity can be given aswhere βTM is another fit parameter. To calculate terminal velocities for variations in fluid viscosity and number density and sphere radius and temperature, a reference state is required where uTM is known for given values of μ, ρ#, R, and T, so the fit parameter, βTM*, can be determined. A reference state is also required to determine the fit parameter, βT*, related to eq . Terminal velocities can then be predicted for variations in the same parameters asTo compare predicted terminal velocities vs particle size from eq to experimental data, results reported by Gibbs et al.[9] were used and are plotted in Figure .

Figure 4

Calculated terminal velocities using Stokes’ law as the resistance term only, uT, and the modified version, uTM, from eq vs particle size compared to experimental data points reported by Gibbs et al.[9] The fit parameters used were βT* = 417 and βTM* = 23 055 and were identified using the following reference state: diameter of 1670 μm with the corresponding velocity of 25.07 cm/s. The estimated uncertainty in the measured data is approximately 4%.

The reference state for determination of fit parameters was a diameter of 1670 μm with a corresponding velocity of 25.07 cm/s, which gave βT* = 417 and βTM* = 23 055, using other input parameters, as given in Table . The expression for the modified terminal velocity, uTM, gives a straight line vs diameter in the range from 50 to 1670 μm and fits measured terminal velocities closely and significantly better than the velocities predicted using the conventional expression for terminal velocity based on Stokes’ law only. The latter curve has the expected parabolic shape vs radius (uT ∼ R2), which follows from eq . A straight line between terminal velocity and particle diameter cannot occur when using the conventional expression for the terminal velocity. The two measured data points for diameters >1670 μm deviate from the general trend probably because of relatively high terminal velocities, which can induce higher flow resistance caused by turbulence (Re numbers greater than 1000).

Table 1

Properties for Calculating Terminal Velocities from Equation and the Curves Shown in Figure

quantity	value	quantity	value
Boltzmann’s constant, kB (J/K)	1.38 × 10–23	water density, ρ (kg/m3)	1000
coefficient of restitution, ε	1.0	water viscosity at 20 °C (mPa s)	1.00
acceleration due to gravity, g (m/s2)	9.81	mass water, m (kg)	2.99 × 10–26
density glass spheres, ρS (kg/m3)	2755	water number density, ρ# (m–3)	3.34 × 1028

Figure shows the percentage influence of viscous resistance vs particle radius for four different viscosities. The curves were generated by analyzing the relative magnitude of the two resistance terms in the modified terminal velocity expression in eq using properties of water at 20 °C as the base (1 mPa s). The three other cases were obtained by changing the numerical value of the viscosity term to 0.01, 100, and 10 000 mPa s without altering other molecular properties. Figure shows that thermal resistance is dominant for particles with radii larger than nanometer-size for a viscosity of 1 mPa, which is the reason for the linear relationship calculated from the modified terminal velocity expression in Figure . Linear relationships have also been reported by others, e.g., Ali et al.[10] and particularly Rubey.[11] To fit experimental data for a wide range of particle sizes, Rubey had to add an additional resistance term, not dependent on viscosity, in an ad hoc manner to Stokes’ law. Hence, an extension of the total flow resistance with an additional term was proposed almost 100 years ago, dictated by the need to obtain satisfactory correspondence between observed and calculated terminal velocities for a wide range of particle sizes. The theoretical justification for a total flow resistance including two terms is suggested to be represented by eq .

Figure 5

Percentage viscous resistance vs particle radius for relative motion of water to spherical particles plotted for four different viscosities (denominator in uTM in eq ) at 20 °C. The cases with viscosities of 0.01, 100, and 10 000 mPa s were implemented by shifting the magnitude of the viscosity term numerically keeping all other molecular properties equal to those of water at 20 °C (Table ).

Viscous resistance becomes more pronounced as the particle size becomes smaller and is dominant for particles of molecular size (≤nm) moving through water with a viscosity of 1 mPa s. It also becomes more influential with increasing fluid viscosity demonstrated by the curve for the highest viscosity value. The thermal part is, however, still dominant for objects with radius in the order of 10–5 m, even when moving through “water” with a viscosity of 10 000 mPa s.

Diffusion in Liquids: Modification of the Stokes–Einstein Relationship

There are many reported studies in the literature stating that it is challenging to fit the Stokes–Einstein relationship to measured diffusion coefficients for molecules or ions through solvents, particularly for varying temperatures.[12−16] Einstein’s expression for the diffusion coefficient, which can also be considered as an early version of the FDT,[30] is given as[37]where the mobility term for the particles is assumed to follow Stokes’ law (MS = ζS–1 = 1/6πμR), and βSE is an adjustable parameter not depending on temperature. Its magnitude depends on the argument stated previously related to the validity for using the thermal resistance expression for liquids and the extent of hydrophilic–hydrophobic properties between the diffusing ions/molecules relative to the solvent molecules.[38,39] Replacing the mobility term in eq using the expression in eq for the friction coefficients gives a modified expression, DM, for the diffusion coefficientwhere βSEM is another fit parameter valid under the same conditions as βSE. To compare eqs and 17 to the experimental data, diffusion and self-diffusion coefficients were calculated vs temperature in the range of 298–423 and 273.17–373 K, respectively, and compared to measured data for diffusion of N2 and CO2 in water and self-diffusion of water reported by Cadogan et al.[17] and Easteal et al.[18] The shapes of N2 and CO2 molecules are assumed to be closely spherical, where deviations from sphericity are accounted for using the hydrodynamic radii as the particle sizes.[40] The observed data points for diffusion of N2 and CO2 were calculated as an average of the reported diffusion coefficients at each temperature neglecting the impact of pressures (incompressible water phase assumed), and the results are depicted in Figure .

Figure 6

Calculated diffusion coefficients for N2 (left) and CO2 (middle) from eqs and 17 compared to experimental data points reported by Cadogan et al.[17] Reproduced from Cadogan, S. P.; Maitland, G. C.; Martin Trusler, J. P. Diffusion Coefficients of CO2 and N2 in Water at Temperatures between 298.15 K and 423.15 K at pressures up to 45 MPa. J. Chem. Eng. Data2014, 59, 519. Copyright 2014 American Chemical Society. To the right, the calculated self-diffusion coefficients for water compared to experimental data points reported by Easteal et al.[18] Other input data are given in Tables –3. The estimated uncertainty in the measured data is approximately 2%.

Table 3

Measured Water Viscosity vs Temperature[41]

water temperature, T (K)	water viscosity, μ (mPa s)
273.01	1.7914
283	1.536
293	1.0016
298	0.890
323	0.547
348	0.3774
373	0.2814
398	0.2222
423	0.1824

The fit parameters identified for N2 and CO2 were βSEM = 1.566 and βSE = 1.47 and βSEM = 1.587 and βSE = 1.504, respectively. For self-diffusion of water, βSEM = 1.3665 and βSE = 1.34 were obtained. The fit parameter values were determined so that the calculated diffusion/self-diffusion coefficients coincided with the measured value at the lowest temperature. Other input data are given in Tables –3. It should be noted that the theoretical expression for viscosity contains the absolute temperature.[2,3,5] It is, however, challenging to calculate viscosity values from the first principle due to mathematical complexities related to realistic molecular pair potentials.[2,3] It is therefore common practice, also used herein, to use measured rather than calculated values for water viscosity vs temperature.

The modified diffusion/self-diffusion coefficient, DM, given in eq , closely follows measured data whereas the conventional diffusion/self-diffusion coefficients calculated from Stokes–Einstein (eq ) are overpredicted at higher temperatures. Since the hydrodynamic radii of N2, CO2, and H2O are small (Table ) and comparable to the radius of the water molecules, viscous resistance is expected to contribute significantly to the total flow resistance according to Figure , particularly at the lowest temperatures. This is confirmed as the ratios of viscous and thermal resistances for N2 and CO2 at 298 K are 16.5 and 18.2, respectively, and decrease to 2.8 and 3.1 at 423 K, respectively. The corresponding ratios for self-diffusion of water are 22.6 and 6.4 at 273.17 and 373.16 K, respectively. Hence, the calculated overprediction of the diffusion and self-diffusion coefficient values from the Stokes–Einstein relationship at higher temperatures is caused by the absence of the thermal resistance term.

Table 2

Properties Used to Calculate the Diffusion/Self-Diffusion Coefficients in Figure a

quantity	value
hydrodynamic radius N2 (pm)	182
hydrodynamic radius CO2 (pm)	165
hydrodynamic radius H2O (pm)	132.5

All hydrodynamic radii are from Ismail et al.[40]

Flow in Porous Media: Modification of Darcy’s Law

The empirical Darcy’s law[1,31] for the volumetric single-phase fluid rate, q, through a porous medium is given bywhere K is the absolute permeability, A is the cross-sectional area, ΔP is the differential pressure, and L is the core sample length. It can be derived and justified theoretically from the N–S equation as demonstrated[1] by Whitaker[42] and Neuman.[43] Hence, it can be deduced based on the results herein that the absolute permeability should depend on temperature, which has also been demonstrated both on the laboratory[20−23] and field scales in geothermal wells.[24−26] Researchers have previously derived an expression for the absolute permeability accounting for temperature variations.[19] Since naturally porous media are complex, the exact relationship between the two resistance terms in eq , valid for an isolated spherical particle, cannot be expected to hold. The grains in such media are not spherical, and the complex morphology is expected to induce a velocity field, which increases the influence of viscous resistance compared to a single sphere; cf. Figure . Two fit parameters, W and ΔPT, were therefore introduced, where W is a weight factor between the viscous and the thermal resistance terms and ΔPT describes “the efficiency” of the thermal resistance. Absolute permeability vs temperature could then be calculated from a modified version of Darcy’s law given bywhere K* is a reference permeability value.[19]Equation was first used to determine the magnitude of W and ΔPT at the lowest temperature corresponding to the highest permeability value K* measured in each series. KCal could then be calculated for increasing temperatures and compared to measured permeability values reported by Aruna[20] and Weinbrandt et al.[21] (Table ). The left panel of Figure shows that calculated permeabilities closely follow measured data in the range of 25.6–148.9 °C for Massillon sandstone data. The corresponding permeability values were 523 mD at 25.6 °C, which was reduced to 219 mD at 148.9 °C, equivalent to a reduction of 58%. The reduction in permeability is even more significant for the Boise sandstone material shown in the right panel of Figure . The initial permeability of 1337 mD at 23.9 °C was reduced to 428 mD at 157.2 °C, equivalent to a reduction of 68%. The horizontal dashed lines in Figure represent permeability vs temperature according to Darcy’s law, i.e., no impact of temperature on the ability of the medium to transmit fluids. Hence, the fits between calculated and measured absolute permeabilities for the two data sets analyzed corroborate the inclusion of a temperature-dependent term in the absolute permeability expression given in eq .

Figure 7

Left: measured and calculated absolute permeabilities as a function of decreasing temperature in the range of 148.9–25.6 °C for Massillon sandstone reported by Aruna.[20] The fit parameters identified were W = 0.75 and ΔPT = 6.2 × 10–5 kg/(m s K1/2). Right: measured and calculated absolute permeabilities for increasing temperature in the range of 23.9–157.2 °C for Boise sandstone cores reported by Weinbrandt et al.[21] The fit parameters identified were W = 0.79 and ΔPT = 6.0 × 10–5 kg/(m s K1/2). Reproduced from Standnes, D. C. Implications of Molecular Thermal Fluctuations on Fluid Flow in Porous Media and Its Relevance to Absolute Permeability. Energy Fuels2018, 32, 8024. Copyright 2018 American Chemical Society. The estimated uncertainty in the measured data is approximately 5%.

There are many observations indicating that absolute permeability varies with temperature.[19−26] The issue is, however, challenging since several other effects may occur simultaneously upon heating naturally porous media, which could impact the measured permeability in several ways. Rosenbrand et al.[22] summarized the effects assumed to be the most important on the laboratory scale as follows: (a) thermal expansion, (b) increased compressibility, (c) mineral dissolution/precipitation, (d) changes of the electrical double layer indirectly causing changes in effective porosity, and (e) particle mobilization due to changes in the surface charge of the minerals, which could cause permeability reduction by plugging of the particles in downstream pore throats.

On larger scales, the so-called thermal fracturing of the rock is normally assumed to be the most important effect to account for the improved injectivity observed when injecting cold water into hot formations conventionally performed for geothermal wells.[24−26] If the injected water is hotter than the formation, the injectivity may decrease upon heating the formation. On the contrary, if the formation cools due to cold water injection, then the injectivity usually increases. Siega et al.[25] investigated systematically the impact of the difference in temperature between the formation and the near-well zone on the injectivity index, II. The temperature variation in the near-well zone was changed by injecting water with lower temperatures than those of the formations. The injectivity index was defined aswhere Q is the mass flow rate, PH is the hydrostatic pressure, WHP is the well head pressure, PF is the pressure due to frictional loss, and PFZ is the pressure at the permeable zone. The pressure difference between the pressure in the well at the permeable zone and in the formation far away from the well expresses the permeability of the medium under prevailing conditions. It is assumed that the formation permeability is the variable that is mostly impacted by changes in temperature of all of the terms in eq . The temperature of the near-well zone due to injection of cold water was measured or estimated and subtracted from the original formation temperature, providing the increase in the temperature on the x-axis. The percentage decrease in II is plotted against this increase in Figure for 21 well observations plus corresponding data from three wells from Iceland. The plot shows that II interpreted as variation in absolute permeability consistently decreases with increasing difference in temperature between the near-well zone and the formation. Hence, the ability to inject water is a strong function of the absolute temperature in the formation. The higher the formation temperature, the lower the ability to inject water. A trend line was identified that captures the trend in the observed data with a correlation coefficient of R2 = 0.82. Variation in II interpreted as variation in absolute permeability vs temperature is also plotted as a dotted line calculated from eq . The input parameters used were ΔPT = 6.0 × 10–5 kg/(m s K1/2) and W = 0.45, which, owing to a lack of detailed information about the rock properties, were obtained using a similar value for ΔPT, as shown in Figure , but by adjusting the value of W. The calculated curve follows the trend line and captures the variations in the observed data in a similar manner. Hence, the observed change in II can be understood in terms of the results obtained herein that absolute permeability is a strong function of absolute temperature.

Figure 8

Percentage reduction in injectivity interpreted as variation in absolute permeability vs difference in temperature between the near-well zone and formation temperature replotted from Siega et al.[25] The dotted line is calculated from eq using input parameters ΔPT = 6.0 × 10–5 kg/(m s K1/2) and W = 0.45. The line is a best-fit trend curve for all observed data given by y = −0.0031x2 + 1.0218x with an R2 = 0.82.

This indicates that thermal resistance is also significant on a large scale for the flow of water in the geothermal wells included in the study. It should also be noticed that the change in II is typically reversible if the temperature of the injected water increases for a period and then decreases again. In this case, a decline in II is observed when the injected water temperature increases, which is recovered when the temperature of the injected water decreases again.[24−26] It should be noticed that improvements in II at lower temperatures are observed even though the viscosity of water increases significantly when the temperature decreases, more than a factor of 6 between 0 and 100 °C. Darcy’s law predicts an inverse relationship between II and viscosity. Hence, the unexpected injectivity behavior of geothermal wells upon alternating injection of cold and hot water can therefore be understood and predicted in terms of the thermal resistance term without the need for additional hypotheses regarding fracture opening/closure due to thermal rock mechanical effects.

d’Alembert’s Paradox

d’Alembert’s paradox refers to the calculated prediction of zero drag on an object located in an incompressible fluid stream, resulting from Euler equations[1,3] or the dimensionless simplified version of the N–S equation in the limit of vanishing fluid viscosity.[32] It is a paradox since experimental data show that a finite drag force acts on the object in the limit of vanishing viscosity. The prediction was calculated by d’Alembert in 1752 based on the potential theory for inviscid fluids. The general solution for the thermal resistance term in one dimension given in eq was approximated by eq in the limit of low relative velocities. The corresponding flow resistance in the limit of high relative velocity, u1 ≫ VT, can also be derived from eq and is given byThe resistance only depends on the mass and the number density of the fluid molecules and the square of the relative velocity. The thermal resistance given by eq for a rectangular object moving through air at 20 °C in one dimension is plotted against relative velocity in Figure using average air property values; therefore, VT = 290 m/s. It can be observed that the relationship between the drag force per square meter and relative velocity is linear below VT (eq ) and proportional to the square of the relative velocity for velocities above VT (eq ). No effects related to turbulent flow behavior are accounted for in Figure . Based on these results, it is therefore proposed that the thermal resistance term alternatively can explain and resolve d’Alembert’s paradox.[32] It predicts a finite resistance to relative motion between fluids and objects both at low and high relative velocities in the limit of vanishing fluid viscosity, as illustrated in Figure .

Figure 9

Drag force per square meter due to thermal resistance as a function of relative velocity between air and a square-shaped object given by the analytical expression in eq .[33] The dotted line shows the magnitude of VT = 290 m/s with a corresponding drag force of 373 kN/m2. Data used are A⊥ = 1 m2 and ε = 1. Fluid data for air are T = 293 K, molecular weight = 29 g/mol, and density = 1.2 kg/m3. No effects related to turbulence are accounted for in the graph.

Discussion

The most important consequence related to the thermal resistance is perhaps that the conventional viscous resistance is negligible for macroscopic objects moving slowly relative to either air or water at around 20 °C. Under such conditions, the thermal resistance is of the order 105 and 106 larger than the viscous resistance for a spherical particle with a radius of 1 cm moving relative to air and water, respectively. Viscous resistance only becomes important when the size of the object approaches molecular dimensions or if the viscosity of the fluids becomes very high. A fluid with waterlike properties must, e.g., have a viscosity of 3 × 106 mPa s to match the magnitude of the thermal resistance term for a spherical object with a radius of 1 cm. Viscous resistance also becomes more important when fluids flow through porous media.[17] It is hypothesized that the increased influence of viscous resistance is caused by the complex morphology and extremely irregular patterns in such media.

Results presented show that modified versions for the terminal velocity of spherical particles falling through liquids, the Stokes–Einstein relationship, and Darcy’s law, all having resistance terms derived from the N–S equation, give improved fit to experimental data compared to the conventional expressions. These observations support the hypothesis of utilizing the “applied approach” stating that eq developed for low-density gases also applies for liquids. They, together with the alternative explanation for d’Alembert’s paradox, corroborate the presence of a thermal resistance term. Introduction of the thermal resistance term is hence able to connect, explain, and predict several observations and trends that cannot be accounted for using viscous resistance as the only flow resistance term. More testing and comparison with measured data are, however, encouraged to further corroborate or falsify the presence of the additional resistance term in wider parameter ranges. This also includes testing of other laws and expressions where fluids and objects are in relative motion, which is not included in the current study.

Conclusions

This work demonstrates that an additional resistance term should be included in the Navier–Stokes equation when fluids and objects are in relative motion. It is based on an observation that the effect of the microscopic molecular random velocity component parallel to the macroscopic flow direction is neglected. The two components of the random velocity perpendicular to the local mean flow direction are accounted for by the viscous resistance, e.g., by Stokes’ law for spherical objects. The relationship between the mean- and the random velocity in the longitudinal direction induces differences in molecular collision velocities and collision frequency rates on the up- and downstream surface areas of the object. This asymmetry therefore induces flow resistance and energy dissipation. The flow resistance resulting from the longitudinal mode is referred to as thermal resistance and depends on the relative velocity between the fluid and the object, the number density and molecular fluctuation statistics of the fluid, and the area of the object and the square root of the absolute temperature. Results show that thermal resistance is dominant compared to viscous resistance considering water and air in slow relative motion to spherical objects larger than nanometer-size at ambient temperature and pressure conditions. Including thermal resistance in the conventional expressions for terminal velocity of spherical particles falling through liquids, the Stokes–Einstein relationship and Darcy’s law, corroborates its presence, as modified versions of these laws fit observed data much more closely than the conventional expressions. The thermal resistance term can alternatively resolve d’Alembert’s paradox as a finite flow resistance is predicted at both high and low relative fluid–object motions in the limit of vanishing fluid viscosity. More testing and reanalysis of the existing experimental data are encouraged to further corroborate or falsify inclusion of the additional resistance term in wider parameter ranges. This also includes testing of other laws and expressions where fluids and objects are in relative motion, which is not included in the current study.