6.7.6.1 Gibbs Free Energy for Electrochemical Reactions

The Gibbs free energy is, of course, also defined for electrochemical reactions where it defines the electromotive force of the reaction

$G=zF\Delta \mu$

If all reagents are provided at STP conditions, ∆G⁰ is given as a function of the standard reduction potential E⁰

$G=zF\Delta E^0$

6.7.6.2 Free Standard Enthalpy of Formation

As the Gibbs free energy depends on the enthalpy, it also suffers from the fact that its values cannot be given in absolute quantity, but only relative to a reference value. As stated, the standard reaction enthalpy ∆H⁰ is referenced to 1 mol of substance at STP conditions while setting ∆H⁰ = 0 for the elements. This gives rise to the free standard enthalpy of formation $\Delta H_f^0$. In an analogy, the Gibbs free energy G of the elements is set G = 0 and the reaction considered at STP conditions for 1 mol of product. This gives rise to the free standard enthalpy of formation

$\Delta G_f^0={\displaystyle \sum \Delta G_{f\text{products}}^0-{\displaystyle \sum \Delta G_{f\text{educts}}^0}}$

Some of the most commonly used $\Delta G_f^0$ values are given in Tab. 6.6.

Tab. 6.6

Free standard enthalpy of formation $\Delta G_f^0$ of important compounds

Compound	$\Delta G_f^0$ kJ mol−1	Compound	$\Delta G_f^0$ kJ mol−1
gases
O3	+163.3	H2O, gas	-228.7
CO	-137.2	CO2	-394.6
SO2	-300.4	NO	86.6
liquids
H2O	70.0	HF	-273.4
HCl	-95.4	HBr	-53.5
HI	1.7	H2O, liquid	-237.3
NH3	-16.5
solids
NaCl	-384.3	CaO	-604.6
Al2O3	-1 583.5	SiO2	-857.3
CuO	-129.8
dissociation
H	+203.4	O	+231.9
F	+62.0	Cl	+105.8
N	+455.9

Free Standard Enthalpy of Formation of Some Reactions. With the values given in Tab. 6.6 the $\Delta G_f^0$ values of chemical reactions can be calculated. As an example, let us consider the following reaction

$\begin{array}{ll}\text{CO}+\frac{1}{2}\text{O}\rightleftharpoons {\text{CO}}_2\Delta G_f^0\hfill & =\Delta G_f^0\left({\text{CO}}_2\right)-\left(\frac{1}{2}\Delta G_f^0\left({\text{O}}_2\right)+\Delta G_f^0\left(\text{CO}\right)\right)\hfill \\ \hfill & =\left(-394.6-\left(0-137.2\right)\right)\text{kJ}{\text{mol}}^{-1}\hfill \\ \hfill & =-257.4\text{kJ}{\text{mol}}^{-1}\hfill \end{array}$

Likewise the reaction of carbon and oxygen gas to carbon dioxide is given by

$\begin{array}{ll}\text{C}+{\text{O}}_2\rightleftharpoons {\text{CO}}_2\Delta G_f^0\hfill & =\Delta {\text{G}}_f^0({\text{CO}}_2)-(\Delta G_f^0(\text{C})+\Delta G_f^0({\text{O}}_2))\hfill \\ \hfill & =-394.6\text{kJ}{\text{mol}}^{-1}\hfill \end{array}$

The reaction of hydrogen gas and chlorine gas is given by

$\begin{array}{ll}\frac{1}{2}{\text{H}}_2+\frac{1}{2}{\text{Cl}}_2\rightleftharpoons \text{HCl}\Delta G_f^0\hfill & =\Delta G_f^0\left(\text{HCl}\right)-\left(\frac{1}{2}\Delta G_f^0\left({\text{H}}_2\right)+\frac{1}{2}\Delta G_f^0\left({\text{Cl}}_2\right)\right)\hfill \\ \hfill & =-95.4\text{kJ}{\text{mol}}^{-1}\hfill \end{array}$

A reaction with a positive value for $\Delta G_f^0$ is the reaction of nitrogen gas and oxygen gas to nitric oxide

$\begin{array}{ll}\frac{1}{2}{\text{N}}_2+\frac{1}{2}{\text{O}}_2\rightleftharpoons \text{NO}\Delta G_f^0\hfill & =\Delta G_f^0\left(\text{NO}\right)-\left(\frac{1}{2}\Delta G_f^0\left({\text{N}}_2\right)+\frac{1}{2}\Delta G_f^0\left({\text{O}}_2\right)\right)\hfill \\ \hfill & =+86.6\text{kJ}{\text{mol}}^{-1}\hfill \end{array}$

6.7.6.3 Interpretation of the Gibbs Free Energy

The Gibbs free energy is an important means for predicting chemical reactions. As already stated, chemical reactions are either exothermic or endothermic, depending on whether or not they require energy from the outside in order to proceed. In adiabatic closed and contained systems, only exothermic reactions can occur as no heat can be transferred from the outside. Transferring heat into a system (endothermic reactions) reduces the entropy of the ambiance while increasing the entropy of the system. Transferring heat out of the system (exothermic reactions) reduces the entropy of the system while increasing the entropy of the ambiance. Therefore, looking at Eq. 6.59 the borderline case is the case for which G = 0 and therefore

$\begin{array}{ll}\Delta S_{\text{system}}\hfill & =\frac{\Delta H}{T}\hfill \\ \Delta S_{\text{ambiance}}\hfill & =-\Delta S_{\text{system}}\hfill \end{array}$

For an exothermic reaction ∆H < 0 and therefore the entropy change of the ambiance is positive as heat is transferred out of the system. For endothermic reactions ∆H > 0 and therefore the entropy change of the ambiance is negative as heat is transferred into the system. However, in chemical reactions, the entropy change due to the reaction must also be considered (see Eq. 6.43). Thus adding up the entropy contribution and the enthalpy change, three cases must be considered.

• $\Delta G_f^0<0$: the reaction will occur spontaneously and work may be gained from it

• $\Delta G_f^0>0$: the reaction will not occur spontaneously; if the entropy is increased (by external heating), the Gibbs free energy may become sufficiently negative for the reaction to occur

• $\Delta G_f^0=0$: the reaction is in equilibrium; usually values between − 5 kJ mol⁻¹ and + 5 kJ mol⁻¹ are considered to be equivalent to 0

It is important to note that merely judging a reaction exothermic and endothermic by looking at the enthalpy changes alone is not a sufficient manner in which to judge if a reaction will effectively occur. As the changes in entropy (especially if the formation or consumption of a gas is involved) may be significant, exothermic reactions may still yield a positive Gibbs free energy and thus, will not occur spontaneously. Also, $\Delta G_f^0$ gives no information about the speed of the reaction, but merely about the degree of completion. Reactions with values $\Delta G_f^0\ll 0$ will tend to run to completion while reactions with values $\Delta G_f^0\gg 0$ will yield almost exclusively the educts. $\Delta G_f^0$ therefore only refers to the thermodynamic stability of a product that is equivalent to considering the equilibrium of a reaction. As an example, we have calculated in section 6.7.6.2 that the reaction of hydrogen gas and chlorine to hydrogen chloride has a negative value $\Delta G_f^0$. However, this gas mixture is stable at ambient conditions. So thermodynamically, the reaction can occur (in principle) - however, it is so slow that it cannot be observed effectively. It is therefore important to note that the Gibbs free energy will indicate if a certain chemical reaction is thermodynamically possible. However, it gives no information on the speed of the reaction.

Determining the Gibbs energy of a substance allows characterizing the ability of this substance to further react. This is sometimes referred to as thermodynamic potential. A system with high potential is able to react further, and a system with low potential will most likely not be reactive. However, the fact that a certain state may be thermodynamically more stable than another does not allow us to directly deduce whether or not this transition really occurs and at what transfer rate. A good example for this is the transition between different states of carbon. Thermodynamically, the graphite state of carbon is more stable than the diamond state. However, at ambient conditions, the rate of this reaction is so low that it cannot be observed.

Diamonds Are Forever? Another interesting example of the interpretation of the Gibbs free energy is the conversion of diamond to graphite, both of which are allotropes. The Gibbs free energy in this case is

$\begin{array}{ll}\Delta G_f^0\left({\text{C}}_{\text{diamond}}\to {\text{C}}_{\text{graphite}}\right)\hfill & =\Delta G_f^0\left({\text{C}}_{\text{graphite}}\right)-\Delta G_f^0\left({\text{C}}_{\text{diamond}}\right)-T\left(\Delta S^0\left({\text{C}}_{\text{graphite}}\right)-\Delta S^0\left({\text{C}}_{\text{diamond}}\right)\right)\hfill \\ \hfill & =\left(0\text{kJ}{\text{mol}}^{-1}-1.9\text{kJ}{\text{mol}}^{-1}\right)-298.15\text{K}\left(5.7\text{kJ}{\text{mol}}^{-1}-2.4\text{kJ}{\text{mol}}^{-1}\right)\hfill \\ \hfill & =-2.88\text{kJ}{\text{mol}}^{-1}\hfill \end{array}$

As we can see, this value is negative, but not very much smaller than 0. So in general, this reaction can occur in principle thermodynamically but it is so slow that the conversion can be neglected. The reason why $\Delta G_f^0\left({\text{C}}_{\text{diamond}}\right)=0\text{kJ}{\text{mol}}^{-1}$ and $\Delta H_f^0\left({\text{C}}_{\text{diamond}}\right)>0\text{kJ}{\text{mol}}^{-1}$ is due to the fact that graphite is the more thermodynamically stable allotrope of carbon and therefore chosen to be zero for $\Delta H_f^0\left({\text{C}}_{\text{graphite}}\right)=0$ reference.

Qualitative Temperature Dependency. In a moment we will have a look at the pressure and temperature dependency of the Gibbs free energy. But we can already discuss the influences of ∆H⁰ and S⁰ qualitatively. We will briefly discuss the four possible cases.

• ∆H⁰ > 0, S⁰< 0: the Gibbs free energy is positive in all cases, there will be no reaction even at elevated temperatures

• ∆H⁰ < 0, S⁰ > 0: the Gibbs free energy is always zero; at lower temperatures the enthalpy change will dominate at higher temperatures the entropy change

• ∆H⁰ < 0, S⁰ < 0: the Gibbs free energy is negative at lower temperatures, but will become positive if the temperature is chosen too high

• ∆H⁰ > 0, S⁰ > 0: the Gibbs free energy is positive at lower temperatures dominated by the enthalpy; at higher temperatures these reactions can proceed due to the increase in dominance of the entropy

A common example is phase changes that can be characterized using the Gibbs free energy. Gases have higher entropy than liquids, which makes their Gibbs energy lower than that of solids. However, at lower temperatures, this term quickly declines and the solid state becomes thermodynamically more stable. This is why liquids and gases are solid at lower temperature.

Equilibrium Point of Water. Let us consider this as an example. Water at STP conditions in the gaseous state has $\Delta H_f^0\left(g\right)=-242.0\text{kJ}{\text{mol}}^{-1}$ and S⁰ (g) = 188.8 J mol^-1. In the liquid state, it has $\Delta H_f^0\left(l\right)=-286.0\text{kJ}{\text{mol}}^{-1}$ and S⁰ (l) = 70.0 J mol⁻¹. The values can be taken from Tab. 6.3 and Tab. 6.4. The equilibrium point is reached, if the two Gibbs free energies are equal

$\begin{array}{l}G\left(g\right)=\Delta H_f^0\left(\text{g}\right)-T_b{S}^0\left(\text{g}\right)\hfill \\ =\Delta H_f^0\left(\text{1}\right)-T_b{S}^0\left(\text{1}\right)\hfill \\ =G\left(1\right)\hfill \end{array}$

From this equation results T_b = 370.37 K = 97.22 °C. This is a very good approximation for the actual boiling point of water, which is the condition at which liquid water is in equilibrium with water vapor. In this example we used the enthalpy and entropy values for STP conditions, for a more precise calculation, the values should have been corrected with increasing temperature.

6.7.6.4 Temperature Dependency of the Gibbs Energy

Before discussing the temperature and pressure dependence of the Gibbs free energy, we will introduce a couple of more important thermodynamic equations. The first equation helps us link the first law of thermodynamics to the second law of thermodynamics, which we have done already in section 6.7.4. When considering closed and contained systems, the heat flow $\stackrel{.}{Q}$ can be expressed using Eq. 6.43 as

$dQ=TdS$

Using Eq. 6.32 and Eq. 6.33 we can deduce

$\begin{array}{l}dW+dQ=dU\hfill \\ -pdV+TdS=dU\hfill \\ -pdV+TdS=d(H-pV)\hfill \\ -pdV+TdS=dH-pdV-Vdp\hfill \end{array}$

$dH=TdS+Vdp$

Both Eq. 6.61 and Eq. 6.62 are used often in thermodynamics and are commonly referred to as fundamental equations. We will also employ them during the derivation of the temperature and pressure dependence of the Gibbs free energy. Now we will write Eq. 6.59 as a differential of the thermodynamic state variables which yields

$\begin{array}{l}dG=dH-d(TS)\hfill \\ dH=dG+d(TS)\hfill \end{array}$

Using Eq. 6.63 we obtain

$\begin{array}{ll}TdS+Vdp\hfill & =dG+d(TS)\hfill \\ dG\hfill & =Vdp-SdT\hfill \end{array}$

We are interested in the change of the Gibbs free energy over temperature for which we have to find the partial differential with respect to the temperature. We will assume the pressure to remain constant, i.e., the reaction is carried out in a sufficiently large vessel. We rearrange Eq. 6.59 to read

$-S=\frac{G-H}{T}$

Using Eq. 6.63 we now derive ${\frac{\partial G}{\partial T}|}_p=-S$. We will need this expression in a moment

${\frac{\partial G}{\partial T}|}_p=-S$

This equates to Eq. 6.64 and thereby yields

$\begin{array}{ll}{\frac{\partial G}{\partial T}|}_p\hfill & =-\frac{H}{T}\hfill \\ {\frac{\partial G}{\partial T}|}_p-\frac{G}{T}\hfill & =-\frac{H}{T}\hfill \end{array}$

Now the left-hand side of this equation can also be written as $T{\frac{\partial G}{\partial T}|}_p$. This is not straightforward to see, but can easily be checked by applying Eq. 3.13

$\begin{array}{ll}{T\frac{\partial \frac{G}{T}}{\partial T}|}_p\hfill & =T\left(\frac{1}{T}{\frac{\partial G}{\partial T}|}_p-\frac{G}{T^2}\right)\hfill \\ \hfill & ={\frac{\partial G}{\partial T}|}_p-\frac{G}{T}\hfill \end{array}$

With this we can rewrite Eq. 6.65 as

$\begin{array}{ll}{\frac{\partial \frac{G}{T}}{\partial T}|}_p\hfill & =-\frac{H}{T}\hfill \\ {\frac{\partial \frac{G}{T}}{\partial T}|}_p\hfill & =-\frac{H}{T^2}\hfill \end{array}$

This equation is referred to as the Gibbs-Helmholtz equation which gives the temperature dependency of the Gibbs energy as a function of the temperature at constant pressure. It was derived independently by Gibbs and German physicist Hermann von Helmholtz.¹ Before solving this differential equation in the following section we will shortly formulate the pressure dependency of the Gibbs free energy.

6.7.6.5 Pressure Dependency of the Gibbs Free Energy

Now we will study the pressure dependency of the Gibbs energy at constant temperature. We employ Eq. 6.63 setting SdT = 0 as the temperature is constant yielding

$\begin{array}{ll}{\frac{dG}{dp}|}_T\hfill & =Vdp\hfill \\ \hfill & =nRT\frac{dp}{p}\hfill \end{array}$

Here we have used the ideal gas equation (Eq. 6.4). It is solved to

$\Delta G_{\text{1}\to \text{2}}=nRT\text{ln}\frac{p_2}{p_1}$

This equation is used in order to calculate the Gibbs energy if the reagents are provided at different partial gas pressures, as the increase in Gibbs energy can simply be calculated as a thermodynamic process 1 → 2. So if the Gibbs energy at STP conditions, i.e., pressure p₀ is given, the Gibbs energy at any pressure can be calculated as

$\begin{array}{l}{G}_p=\Delta G^0+\Delta G_p{}_{0\to p}\hfill \\ =\Delta G^0+nRT\ln \frac{p}{p_0}\hfill \end{array}$

If more than one reaction partner takes part in the reaction, the changes in Gibbs energy are summarized. If we now consider the arbitrary reaction

${\text{n}}_{\text{A}}\text{A }+{\text{ n}}_{\text{B}}\text{B}\rightleftharpoons n_{C}C$

where the coefficients n_i account for the stoichiometry of the reaction, the overall change of Gibbs free energy can be obtained as the sum of the changes of the individual compounds as

$\begin{array}{ll}\Delta G_{\text{start}\to \text{end}}\hfill & =\Delta G^0+RT\left(n_C\ln \left(\frac{p_{C,\text{end}}}{p_{C,\text{start}}}\right)-\left(n_A\ln \left(\frac{p_{A,\text{end}}}{p_{A,\text{start}}}\right)+n_B\ln \left(\frac{p_{B,\text{end}}}{p_{B,\text{start}}}\right)\right)\right)\hfill \\ \hfill & =\Delta G^0+RT\left(\ln \left(\frac{p_{C,\text{end}}^{n_C}{p}_{A,\text{start}}^{n_A}{p}_{B,\text{start}}^{n_B}}{p_{C,\text{start}}^{n_C}{p}_{A,\text{end}}^{n_A}{p}_{B,\text{end}}^{n_B}}\right)\right)\hfill \\ \hfill & =\Delta G^0+RT\ln \left(\frac{K_p}{K_p^0}\right)\hfill \end{array}$

where we have introduced the equilibrium constant of the reaction at a given pressure K_p and the equilibrium constant of the reaction at STP conditions $K_p^0$. The equilibrium can also be described for solutions using the respective molar concentrations at STP conditions K_c and $K_p^{\mathrm{c}}$.

If the reaction is carried out in equilibrium and in a very large vessel, the Gibbs free energy does not change and ∆Gstart→end = 0. In this case, ∆G⁰ can be calculated as

$\Delta G^0=-RT\ln \left(\frac{K_p}{K_p^0}\right)$

Electrochemical Reactions. For electrochemical reactions, the dependency of the Gibbs energy on the concentrations is likewise described by

$\Delta G_{\text{start}\to \text{end}}=E^0-\frac{RT}{zF}\ln \left(\frac{K_p}{K_c^0}\right)$

Van’t Hoff Equation. We now finish the formulation of the temperature dependency of the Gibbs free energy. We have already derived the Gibbs-Helmholtz equation (see Eq. 6.66) which we rewrite to

$\frac{d}{dT}\frac{\left(\Delta G^0\right)}{T}=-\frac{\Delta H^0}{T^2}$

where we use Eq. 6.68 for the left-hand side to obtain

$\frac{d}{dT}\ln \left(\frac{K_p}{K_p^0}\right)=\frac{\Delta H^0}{R{T}^2}$

This equation is referred to as the Van’t Hoff¹equation. It is an important equation in practical chemistry as it forms the basis of the Q₁₀temperature coefficient which predicts that the reaction speed of a chemical reaction is approximately doubled (sometimes even tripled) if the temperature is raised by 10 K. This prediction is surprisingly correct for many systems in chemistry, biochemistry, and biology.

Equation Eq. 6.69 can be solved using Eq. 6.68 to result in

$\ln \left(\frac{K_p}{K_p^0}\right)=\frac{\Delta H^0}{R}\left(\frac{1}{T}-\frac{1}{T_0}\right)$

This equation gives the dependency of the reaction equilibrium with respect to temperature. If the equilibrium constant K_p1 at a temperature T₁ is known, the equilibrium constant K_p2 at a temperature T₂ can be calculated.

We have thereby established an analytical correlation between the Gibbs energy and the pressure as well as the temperature, respectively. As stated, the Gibbs energy is dependent on these parameters which we may change during the reaction. If changing these parameters would result in a shift of the Gibbs free energy, this may help optimizing reaction conditions.

6.9.2 Fourier’s Law of Heat Conduction

6.9.2.1 Stationary Heat Conduction

Fourier’s Law of Heat Conduction in One Dimension. Fourier’s law of heat conduction is named after French mathematician Jean-Baptiste Fourier. The history of this important fundamental equation of heat conduction has been reviewed by Narasimhan [14] in an article worth reading. Fourier observed that the heat transported through a thin layer of a material was proportional to the differences in temperature of the two sides of the layer (Fig. 6.5a)

$\begin{array}{l}{\dot{q}}_A=-k\frac{T_2-T_1}{X}\hfill \\ =-\lambda \frac{T_2-T_1}{X}\hfill \\ =-\lambda \frac{dT}{dx}\hfill \end{array}$

As the specific heat ${\dot{q}}_A$ is proportional to the area across which the heat is transferred, it makes sense to refer the heat to the area therefore

$d\dot{Q}={\dot{q}}_{A}dA$

$=-\lambda \frac{dT}{dx}dA$

The unit of ${\dot{q}}_A$ is therefore J m⁻² s⁻¹. Eq. 6.71 is Fourier’s law of heat conduction in one dimension. The negative sign indicates that the heat flux is directed from the warmer to the colder end of the layer. The proportionality factor is referred to as the material’s thermal conductivity λ. This is a material constant and independent of the geometry. However, it is dependent on the temperature although in most fluid mechanical calculations, it is often assumed to be a constant. Tab. 6.8 lists a selection of thermal conductivity values for solids. Tab. 9.3 lists a selection of thermal conductivity values for liquids and gases.

Fourier’s Law of Heat Conduction in Three Dimensions. The temperature profile that is formed due to this difference in temperature must be calculated. It can be found by extending Eq. 6.71 to the infinitesimal volume dV (Fig. 6.5b). We will derive Fourier’s law of heat conduction in three dimensions. The temperature changes along the axes are given by

in x-direction

$\begin{array}{l}d{T}_x=T_{x+\Delta x}-T_x\hfill \\ =T_x+\frac{\partial T}{\partial x}dx-T_x\hfill \\ =\frac{\partial T}{\partial x}dx\hfill \end{array}$

in y-direction

$\begin{array}{l}d{T}_y=T_{y+\Delta x}-T_y\hfill \\ =T_y+\frac{\partial T}{\partial y}dy-T_y\hfill \\ =\frac{\partial T}{\partial y}dy\hfill \end{array}$

in z-direction

$\begin{array}{l}d{T}_z=T_{z+\Delta x}-T_z\hfill \\ =T_z+\frac{\partial T}{\partial z}dz-T_z\hfill \\ =\frac{\partial T}{\partial z}dz\hfill \end{array}$

Substituting Eq. 6.74, Eq. 6.75, and Eq. 6.76 into Eq. 6.71 allows us to rewrite

$\begin{array}{ll}{\dot{q}}_{A,\text{total}}\hfill & =-\left(\lambda \frac{\partial T}{\partial x}dx+\lambda \frac{\partial T}{\partial y}dy+\lambda \frac{\partial T}{\partial z}dz\right)\hfill \\ \hfill & =-\lambda \left(\overrightarrow{\nabla }\cdot T\right)\hfill \end{array}$

that is the Fourier law for heat conduction is three dimensions. The operator $\overrightarrow{\nabla }$ is the nabla operator which will be introduced in section 7.1.3.1.

Balance of Heat Flows. We will now turn to the balance of heat flows in the infinitesimal volume (Fig. 6.5c). It is given by

in x-direction

$\begin{array}{ll}d{\dot{q}}_{A,x}dydx\hfill & ={\dot{q}}_x-\left({\dot{q}}_x+\frac{\partial {\dot{q}}_x}{\partial x}dx\right)dydz\hfill \\ \hfill & =-\frac{\partial {\dot{q}}_x}{\partial x}dxdydz\hfill \end{array}$

in y-direction

$\begin{array}{ll}d{\dot{q}}_{A,y}dxdy\hfill & ={\dot{q}}_y-\left({\dot{q}}_y+\frac{\partial {\dot{q}}_y}{\partial y}dy\right)dxdz\hfill \\ \hfill & =-\frac{\partial {\dot{q}}_y}{\partial y}dxdydz\hfill \end{array}$

in z-direction

$\begin{array}{ll}d{\dot{q}}_{A,z}dxdy\hfill & ={\dot{q}}_z-\left({\dot{q}}_z+\frac{\partial {\dot{q}}_z}{\partial z}dz\right)dxdy\hfill \\ \hfill & =-\frac{\partial {\dot{q}}_z}{\partial z}dxdydz\hfill \end{array}$

in summary

$\begin{array}{ll}d{e}_{V,\text{tot}}dxdydz\hfill & =-\left(\frac{\partial {\dot{q}}_x}{\partial x}+\frac{\partial {\dot{q}}_y}{\partial y}+\frac{\partial {\dot{q}}_z}{\partial z}\right)dxdydz\hfill \\ d{e}_{V,\text{tot}}\hfill & =-\left(\frac{\partial {\dot{q}}_x}{\partial x}+\frac{\partial {\dot{q}}_y}{\partial y}+\frac{\partial {\dot{q}}_z}{\partial z}\right)\hfill \end{array}$

where we have used the first law of thermodynamics (see Eq. 6.32) in order to relate the heat flows to the change of energy. Using Eq. 6.71, Eq. 6.78 can be rewritten as

$d{e}_{V,\text{tot}}=\lambda \left(\frac{{\partial }^{2}T}{\partial x^2}+\frac{{\partial }^{2}T}{\partial y^2}+\frac{{\partial }^{2}T}{\partial z^2}\right)$

$=\lambda \Delta T$

Here we have used the Laplace operator ∆ which will be introduced in section 7.1.3.5.

Stationary Heat Conduction.Stationary heat conduction refers to the fact that we do not take into account warming of the fluid, i.e., there is no change of enthalpy. If the control volume does not have any additional heat sources (e.g., en exothermic chemical reaction) or sink (e.g., a endothermic reaction) the change of total energy Eq. 6.80 becomes a classical Laplace equation (see section 8.1.6)

$\lambda \Delta T=0$

6.9.2.2 Instationary Heat Conduction

Instationary heat conduction takes into account changes of enthalpy of the fluid control volume. We therefore have to use Eq. 6.80 and substitute the left-hand side of the equation by the change of enthalpy over time as a consequence of heat conduction. This equation results in

$\begin{array}{ll}\rho c_p\frac{\partial T}{\partial t}\hfill & =\lambda \Delta T\hfill \\ \frac{\partial T}{\partial t}\hfill & =\frac{\lambda }{\rho c_p}\Delta T\hfill \\ \frac{\partial T}{\partial t}\hfill & =\alpha \Delta T\hfill \end{array}$

where we have introduced the thermal diffusivity α which has the unit mm² s⁻¹. It is composed only of material constants and therefore is also a constant for a given material (see section 9.9.3). Tab. 9.4 lists a selection of thermal diffusivity values for commonly encountered liquids. Eq. 6.82 is the Fourier law for instationary heat conduction in three dimensions.

6.9.3 Diffusion

6.9.3.1 Introduction

Diffusion is the random movement of small particles over time. This effect of the movement of small particles across large distances was first described analytically by German physician Adolf Fick.¹ It is one of the two important mass transport phenomena that we encounter in fluid mechanics (the other being convection). Convection is the term used for the bulk movement of liquids, e.g., in the form of flowing fluids. If the fluid is at rest, no mass will be transported. On the other hand, diffusion is the mass transport phenomena that can occur even when the bulk of the liquid is at rest. It is based on the action of Brownian motion which is sufficient to move particles if they are sufficiently small, or light. Therefore, diffusion is only relevant for molecules of sufficiently small size.

6.9.3.2 Visualizing Diffusion: A Digital Diffusion Experiment

Rationale. Many students have difficulty with the concept of diffusion because it seems so different from the Newtonian mechanics we are accustomed to. Usually particles move because they have been exposed to forces, thus they gained momentum which essentially conserved the action of the forces that accelerated them. This type of movement seems intuitive: If we apply a force, a particle will move. Diffusion seems to work on entirely different mechanics which is not entirely true.

Diffusion is actually a very simple process which can be visualized very easily. We will do this using a very simple Maple listing which is shown in listing 6.1. This is in the true sense of the meaning a “digital diffusion” experiment, i.e., a computer program that emulates surprisingly exactly what happens during diffusion.

1   restart:read "Core.txt":
2
3   #this function makes the given number of steps using a random number generator
4   dosteps:=proc(generator,start,epsilon,timesteps) local i,traveled;
5    traveled:=0;
6    for i to timesteps do
7     if generator()=1 then
8       #move left
9       traveled:=traveled-epsilon;
10     else
11      #move right
12      traveled:=traveled+epsilon;
13     end if;
14   end do:
15   return start+traveled;
16  end proc:
17
18  #this function creates our diffusion model
19  diffuse:=proc(moleculecount,epsilon,timesteps) local i,generator,output;
20   #create the random number generator
21   generator:=rand(0..1);
22   #this is the output: x-values (first column) and y-values (second column)
23   output:=Matrix(moleculecount,2);
24   for i to moleculecount do
25     #x-value
26     output[i,1]:=dosteps(generator,0,epsilon,timesteps);
27     #y-value
28     output[i,2]:=dosteps(generator,0,epsilon,timesteps);
29   end do:
30   return output:
31  end proc:
32
33  #perform 'digital diffusion'
34  quickplot([diffuse(10,1,10),diffuse(10,1,50),diffuse(10,1,100)],-300..300,-300..300,"x [a.u.]", "y
       [a.u.]",["10 steps","50 steps","100 steps"],true);
35  quickplot([diffuse(10,5,10),diffuse(10,5,50),diffuse(10,5,100)],-300..300,-300..300,"x [a.u.]", "y
       [a.u.]",["10 steps","50 steps","100 steps"],true);
36  quickplot([diffuse(10,10,10),diffuse(10,10,50),diffuse(10,10,100)],-300..300,-300..300,"x [a.u.]",
       "y [a.u.]",["10 steps","50 steps","100 steps"],true);
37  quickplot([diffuse(10,20,10),diffuse(10,20,50),diffuse(10,20,100)],-300..300,-300..300,"x [a.u.]",
       "y [a.u.]",["10 steps","50 steps","100 steps"],true);

Suppose we have a very very small particle, e.g., an ion. If a particle is very small, its mass will also be very small. Even if not accelerated by forces, this particle will eventually start to move due to Brownian motion. Brownian motion is usually neglected because it is very much smaller than Newtonian movement. However, if the particle is very small, the forces associated with Brownian motion can have a certain effect and the particle will start to move. Obviously this motion will increase as the temperature increases. Please note that this movement is very small, especially compared to Newtonian movement. The particle will move around its original position. The extent of this movement will be influenced by a number of factors. The mass and the size of the particle will obviously play a dominant role. Larger and heavier particles will move less compared to lighter and smaller particles. The density and the viscosity of the surrounding also plays a crucial role. In dense and highly viscous liquids, the particle movement is hindered due to higher friction forces. All of these factors can be summarized in a variable we will call . is the distance that a particle moves from its original position in a given time interval ∆t. The longer this time interval, the longer the traveled distance will be.

For the sake of simplicity, we will only consider one-dimensional particle movement, i.e., the particle only travels along one axis (in this case the x-axis). Let us now assume our particle starts at the beginning of our experiment at the origin x (t = 0) = 0 (see Fig. 6.6). Our particle is now free to move either to the left or to the right side. Under perfect conditions, the probability is 0.5 and therefore identical for both sides. Therefore it is a purely random event if the particle will move to x (t = ∆t) = or to x (t = ∆t) = −. In both cases,

the particle will have moved. Now during the second time interval, the particle is again free to move either to the left or right. Therefore it may end up as x (t = 2∆t) = 2 , x (t = 2∆t) = −2 , or x (t = 2∆t) = 0 where the probability for x = 0 is 0.5, and the probability for x = 2 and x = −2 is 0.25, respectively. The chances are thus higher of finding our particle at the origin after two time steps. Here you can see why diffusion has a probabilistic element to it: The particle’s position is a probability distribution. It is simply more likely to be in the origin than anywhere else but it may also have moved.

Maple Worksheet. If you look at listing 6.1 we can simply emulate this probabilistic movement via a random number generator. This is what the function dosteps does (see line 4). It takes a random number generator generator as arguments as well as a particles original location start. It also requires the information how far a particle may travel. This information is given by the argument epsilon or as we have been using. The last argument timesteps indicates the number of sequential steps we want our particles to move, i.e., the number of time intervals ∆t. The function dosteps iterates the number of time steps. For each time step, it requests a number from the random number generator which will only return 0 or 1. In case of it returning 1, the particle moves to left by decrementing the variable traveled by epsilon. If the random number generator returns 0 the variable traveled is increment by epsilon thus the particle moves toward positive x, i.e., to the right. After the given number of steps, the function returns the sum of the original position start and the traveled distance traveled. Therefore this function effectively emulates a particle which can travel the distance per time step and returns to its location after timesteps number of timesteps assuming its original position start.

The second function we define is diffuse (see line 19). This function takes three parameters: the number of molecules moleculecount, epsilon, and timesteps the latter two of which we already know. The function first creates a random number generator that only creates the values 0 and 1 with equal distribution. In the next step, it creates a matrix with two columns and a number of lines equivalent to the number of molecules. It then iterates and, for each molecule, it calls the function dosteps both for the x-value (first column) and for the y-value (second column). The molecules are “seeded”, i.e., placed at the origin x = y = 0 at the start of the experiment.

The output of this “digital diffusion experiment” is shown in Fig. 6.7. In this experiment a total number of 20 molecules were seeded at the origin. Depending on the number of timesteps, i.e., the duration of the experiment and the “diffusion step width” one can clearly see that the cloud of molecules distributes, i.e., smears. This is exactly what happens in diffusion. The plots are created using the code starting from line 19. You will see that every time you run these commands the plots will look slightly differently. This is due to the fact that diffusion is a random process and the output will look slightly different every time.

Diffusion in More Than One Dimension. It may seem strange that we can simply convert from one-dimensional to two-dimensional diffusion simply by calling the function dosteps for x and y. In fact this is correct because during our assumption of diffusion, we assume the particle to be able to move in one direction. Obviously we would need to correct the lengths at which the particle can travel in order to ensure that the total distance traveled (which is $\sqrt{{\left(\Delta x\right)}^2+{\left(\Delta y\right)}^2}$) is equal to but this does not change the overall picture. We may rewrite the algorithm such that the random number generator yields three values which can be interpreted as

• Indicator if the particle movement in x is in the positive or negative direction (first value)

• Indicator if the particle movement in y is in the positive or negative direction (second value)

• Fraction values that determines at which extent will be contributed to movement along the x-axis or the y-axis

6.9.3.3 Conservation of Mass

In the following, we will derive the two laws of Fick for diffusion. We start with the first law which describes the mass flux into and out of a control volume. Here, we use the mass flow normalized to a surface which is generally referred to as the mass flux J with the unit kg m⁻² s⁻¹. Balancing the mass flux at the unit cell (see Fig. 6.8) we find

$\begin{array}{ll}\frac{\partial m}{\partial t}\hfill & =\left(J_x-\left(J_x+\frac{\partial J_x}{\partial x}\right)\right)dydz\hfill \\ \hfill & +\left(J_y-\left(J_y+\frac{\partial J_y}{\partial y}\right)\right)dxdz\hfill \\ \hfill & +\left(J_z-\left(J_z+\frac{\partial J_z}{\partial z}\right)\right)dxdy\hfill \\ \hfill & =-\left(\frac{\partial J_x}{\partial x}+\frac{\partial J_y}{\partial y}+\frac{\partial J_z}{\partial z}\right)dxdydz\hfill \\ \hfill & =--\left(\frac{\partial J_x}{\partial x}+\frac{\partial J_y}{\partial y}+\frac{\partial J_z}{\partial z}\right)dV=-\overrightarrow{\nabla }\overrightarrow{J}dV\hfill \end{array}$

We divide this equation by dV and obtain the conservation of mass as a function of the mass concentration ρ and the mass flux J

$\frac{\partial \rho }{\partial t}=-\overrightarrow{\nabla }\overrightarrow{J}$

where we have again used the nabla operator which will be formally introduced in section 7.1.3.1.

Derivation Using the Gauss’s Theorem. We can derive the same equation using Gauss’s theorem (see section 7.2.1) which gives us the amount of substance transported across the boundary of an arbitrary control volume is given by

${\displaystyle \oint \overrightarrow{J}\cdot \overrightarrow{n}dA={\displaystyle \int \overrightarrow{\nabla }\overrightarrow{J}dV}}$

which must be equal to the negative change of mass inside of the control volume. Therefore

$\begin{array}{ll}-\frac{\partial m}{\partial t}\hfill & ={\displaystyle \int \overrightarrow{\nabla }\overrightarrow{J}dV}\hfill \\ \frac{\partial \rho }{\partial t}\hfill & =-\overrightarrow{\nabla }\overrightarrow{J}\hfill \end{array}$

6.9.3.4 Fick’s First Law of Diffusion

Diffusion is a process that can be derived from kinetic gas theory. It is simply the relation between the kinetic energy of a particle or a molecule, i.e., its inertia and the friction forces of the surrounding fluid. At very low Reynolds numbers (Re 1, see section 9.9.8 for an introduction to the Reynolds number) a particle will experience a friction proportional to its velocity

$\overrightarrow{v}=\alpha \overrightarrow{F}$

where σ is referred to as mobility. The mobility is inverse-proportional to the size of the molecule, i.e., its diameter d and the viscosity η of the fluid: The higher the diameter and the more viscous the fluid, the stronger the mobility will be hindered. The drag that the surrounding fluid will have on the particle is given by Stoke’s¹drag law which states that

$\begin{array}{ll}{F}_{\text{drag},\text{sphere}}\hfill & =6\pi \eta r\hfill \\ \sigma \hfill & =\frac{1}{6\pi \eta r}\hfill \end{array}$

In general, a chemical potential (which has the unit of energy) is dependent on the local concentration of a compound according to

$\mu ={\mu }_0+k_{B}T\ln \frac{c}{c_0}$

If this potential has a gradient, it will create a force according to

$\overrightarrow{F}=-\overrightarrow{\nabla }\mu =-\frac{k_{B}T}{c}\overrightarrow{\nabla }c$

Using Eq. 6.85, Eq. 6.86, and Eq. 6.87 we find

$\begin{array}{ll}\overrightarrow{v}\hfill & =\sigma \overrightarrow{F}=\frac{1}{6\pi \eta r}\overrightarrow{F}\hfill \\ \hfill & =\frac{k_{B}T}{6\pi \eta rc}\overrightarrow{\nabla }c\hfill \\ \overrightarrow{v}c\hfill & =\frac{k_{B}T}{6\pi \eta r}\overrightarrow{\nabla }c\hfill \\ J\hfill & =\frac{k_{B}T}{6\pi \eta r}\overrightarrow{\nabla }c\hfill \\ J\hfill & =D\overrightarrow{\nabla }c\hfill \end{array}$

where we have introduced the massflux J with the unit kg m⁻² s⁻¹ and the diffusion or diffusivity coefficient D as

$D=\frac{k_{B}T}{6\pi \eta r}$

The mass flux is the amount of substance transported through a given surface area per time. D gives a relation of how well a gradient in concentration can be “turned into” mass movement. For low mobility molecules, e.g., heavy or bulk molecules, as well as for highly viscous environments, these values are small. Eq. 6.88 is referred to as Fick’s first law of diffusion. It states that the mass flux is proportional to the gradient in concentration.

Returning to our experiment with the “digital diffusion”, this is surprisingly straightforward. We assume the particle movement to be equally probable in the positive and the negative x-direction. However, assume there is a larger conglomeration of particles in the negative x-direction. As the volume is more “crowded” in this direction, the particle movement in the other direction is more probable. The movement of the particle is hindered by the higher density of particles in the negative x-direction. An everyday analogy to this is a traffic jam: The more cars, the slower the movement and in consequence, fewer cars pass by per unit time. Obviously, the direction of the particle movement is in the opposite direction of the concentration gradient, i.e., away from it (thus the negative sign).

The proportionality factor for this relation is generally referred to as the diffusion constant D. It has the unit m² s⁻¹ which is chosen such that the units of the equation match up. Referring to our “digital diffusion” experiment, it may be interpreted as a combination of the time step ∆t and the stepwidth . The higher this value, the faster a molecule can move. Obviously, mass transport will be significantly faster at higher diffusion constants.

Diffusion Time. Mass transport by diffusion is not a fast process. The root-mean-square position of a particle after a given time interval t, i.e., the distance the particle has traveled in this time interval is given by

$\langle \overline{x}\rangle =\sqrt{2tD}$

The mathematics of this relation is derived in section 8.3.8.5. In general the distance traveled by molecule scales with the square root of both the time and the diffusion coefficient. Tab. 6.9 lists some commonly encountered diffusion coefficients. Fig. 6.9 displays the difference in mobility of molecules of different size. As you can see, mass transport by diffusion becomes a very slow process for bulkier molecules.

6.9.3.5 Fick’s Second Law of Diffusion

Fick’s second law of diffusion combines the mass conservation equation with the first law of Fick (see Eq. 6.84 and Eq. 6.88). It yields the relationship between the gradient in the mass concentration and the change of mass as

$\begin{array}{ll}\frac{\partial \rho }{\partial t}\hfill & =D{\overrightarrow{\nabla }}^2\rho \hfill \\ \hfill & =D\Delta p\hfill \end{array}$

Again, the Laplace operator ∆ will be formally introduce in section 7.1.3.5. Eq. 6.90 is a homogeneous second-order PDE which is simply referred to as the diffusion equation. In section 8.3.8 we look into methods of solving this equation.

6.9.3.6 Combined Convection and Diffusion

In many practical fluid mechanical problems, diffusion is often neglected due to the slow mass transport over long length scales. Mass transport is therefore solely based on the in- and outflow of mass into and out of the control volumes. This fluid movement is generally referred to as convection. However, in many microfluidic applications, the length scales are reasonably small and it may be necessary to account for both diffusion and convection. We have already derived Fick’s second law of diffusion (see Eq. 6.90) which is the most general PDE for describing diffusion.