Thermodynamic analysis of thermal convection based on entropy production.

Flow patterns have a tendency to break the symmetry of an initial state of a system and form another spatiotemporal pattern when the system is driven far from equilibrium by temperature difference. For an annular channel, the axially symmetric flow becomes unstable beyond a given temperature difference threshold imposed in the system, leading to rotational oscillating waves. Many researchers have investigated this transition via linear stability analysis using the fundamental conservation equations and the generic model amplitude equation, i.e., the complex Ginzburg-Landau equation. Here, we present a quantitative study conducted of the thermal convection transition using thermodynamic analysis based on the maximum entropy production principle. Our analysis results reveal that the fluid system under nonequilibrium maximizes the entropy production induced by the thermodynamic flux in a direction perpendicular to the temperature difference. Further, we show that the thermodynamic flux as well as the entropy production can uniquely specify the thermodynamic states of the entire fluid system and propose an entropy production selection rule that can be used to specify the thermodynamic state of a nonequilibrium system.

Introduction

Evaluating the stability of nonequilibrium states and finding a transition point between two nonequilibrium states require the construction of governing equations and subsequent detailed and laborious analysis of the equations[1-4]. Variational principles are quite significant for the analysis of phase behavior under nonequilibrium processes because of their broad application and the ease with which they are handled. The results of theoretical and experimental studies have suggested that nonequilibrium states maximize the entropy production of nonequilibrium processes (the so-called maximum entropy production principle (MEPP))[5-11]. The variational principles obtained from the MEPP[5] yield well-established equations and various constitutive equations, such as the diffusion equation, Leslie–Ericksen equation, Nernst–Fokker–Planck equation, and constitutive equations describing complex phenomena, such as thin film hydrodynamics, colloid particles, kinetics of phase separation, deformation and diffusion in gels, nonlinear viscoelasticity of polymers, liquid crystals, and drying droplets[12-16].

MEPP-based analysis has been extensively discussed for prediction of a transition point between two nonequilibrium states in complex systems; for example, those found in the configurational changes of crystal growth and the mode changes in droplet oscillation, which involves the two nonequilibrium processes interfering with each other, i.e., mass transfer and heat conduction, and mass transfer and viscous dissipation, respectively[6,7,11]. Nonequilibrium processes can be divided into two types: compound and complex[10]. For compound processes, the uncoupled processes are decomposed into various elementary processes that are only dependent on the corresponding thermodynamic flux, and not on all of the fluxes. All of the entropy production in compound processes can be represented as the sum of the functions of all of the elementary processes. For complex processes, various elementary processes are coupled and interfere with each other. Previous studies that attempted to disprove the MEPP are outside the range of MEPP applicability because the nonequilibrium systems investigated are compound systems[17-22]. It is claimed that the MEPP is valid for complex systems and invalid for compound systems[23]. The MEPP indicates that the systems evolve such as to maximize their entropy production. However, it is unclear whether partial or total entropy production is an essential factor to determine the time course of the systems. Further, previous studies have not explained why a certain component of entropy production can provide a selection rule that determines the nonequilibrium states although different entropy production components originate from various irreversible processes occurring within a nonequilibrium system. The recent literatures on convection experiments demonstrate that the realized convection patterns are not governed by principles like the MEPP[24-26]. Thus, it was stated that the MEPP has been applied in a largely ad hoc manner and its successes remain something of an unexplained curiosity[27].

In this study, we evaluated the relationship between the flow patterns and entropy production via numerical simulation to develop an entropy production selection rule that can be used to specify the thermodynamic state of a nonequilibrium system. We consider the flow patterns driven by thermocapillary instabilities as a complex process involving two irreversible processes: viscous dissipation and heat conduction. The temperature difference imposed across the cell induces a surface tension gradient on the free surface of the fluid, leading to a surface flow towards the cold side. The resulting thermal convection induces a various flow patterns with an internal circulation. In the case of an annular channel, the flow pattern is axially symmetric along the temperature gradient with an internal circulation. This axially symmetric flow (ASF) becomes unstable beyond a given temperature difference threshold and subsequently symmetry-breaking flow, i.e., rotational oscillating waves, appears. The oscillating waves propagate in a direction perpendicular to the temperature gradient applied to the system, i.e., in the circumferential direction, and the temperature changes periodically. This rotational oscillating flow is called hydrothermal wave (HTW)[28-30]. The research of the transition from ASF to HTW has been stimulated by numerical experiments on the fundamental conservation equations, i.e., the Navier-Stokes equation and conservations of mass and energy, and the generic model amplitude equation, i.e., the complex Ginzburg-Landsu equation. Remarkable predictions concern the existence of various flow patterns and instabilities[31,32]. In this study, we first computed each component of thermodynamic flux and thermodynamic force in both ASF and HTW to calculate each component of entropy production. Then, we ascertained which component of thermodynamic flux can specify the system, and which component of entropy production can be maximized to predict the system behavior.

We assume that a linear relationship between the thermodynamic force and the thermodynamic flux holds for the two irreversible processes considered in this study. Then, for heat conduction, the thermodynamic force is XT = ∇(1/T), and the thermodynamic flux is JT = −λ∇T, where λ is the thermal conductivity and T is the temperature; for viscous dissipation of the fluid, they are XV, = σ/T and JV, = ∂u/∂x, respectively, where σ is the viscous stress tensor, u is the i-th component of the fluid velocity vector , and x is the j-th component of coordinates. The entropy production of each elementary process can then be calculated as the product of the thermodynamic force and thermodynamic flux. We find that the ratio of both entropy productions, i.e., for heat conduction (σT) and viscous dissipation (σV), is σV/σT ≈ 10−7. The energy dissipation due to the viscosity is much less than the entropy production due to heat conduction; therefore, we ignore the energy dissipation in this study. Note that this question of ignoring the smaller entropy production requires careful consideration and will be discussed in a later work. Further, in the analysis presented below, we omit the subscripts T and V from the physical quantities because we focus on the entropy production due to heat conduction only.

We take the time-averaged values of the local thermodynamic flux during a certain period τ as and the local thermodynamic force as where i = x, y, z, because the thermodynamic variables vary spatiotemporally in the nonequilibrium system. Most variables are measured at the point (x, y, z) = (30 mm, 0 mm, 3 mm). The x component of the heat flux is in the same direction as ΔT between the inner and outer walls. We calculate the absolute value of the local entropy production from the following relation,where JX denotes the simple product of J and X. Note that the entropy production is scalar, and σ represents the entropy production generated by the ith component of the heat flux. The total entropy production is the sum of the entropy production contributions for all components:

Figure 1 shows the computed surface temperature fluctuation, δT, in the horizontal plane, which is defined by When the temperature difference ΔT is less than 7 K, an axially symmetric flow (ASF) occurs. However, as ΔT increases above 8 K, traveling waves propagate in the direction perpendicular to the temperature gradient applied to the system, i.e., in the circumferential direction, and the temperature changes periodically. This rotational oscillating flow is called an hydrothermal waves (HTW)[28-30]. In a temperature range for which an HTW is produced, the surface temperature started to vibrate at about 50–100 s after the start of the simulation. The amplitude of the temperature oscillation approached a constant value over time, and, after 400 s, the wave number of the transmission wave also became a constant value with the standard deviation of the amplitude being less than 0.01 K. However, even if the simulation was performed for 2500 s, the amplitude of the temperature oscillation did not become constant, and there was a standard deviation of 0.0008 K. Considering the calculation cost, we analyzed the two intervals of 200–300 s and 300–400 s in this study.

Figure 1

Schematics of the annular container (upper panel). Surface temperature difference at various temperature differences and at t = 150 s (lower panels).

Results and Discussion

First, we investigated the relationship between |J| and |X| (Fig. 2a). The relationship is a single line through the origin, although the flow pattern changes from ASF to HTW. Moreover, the relationship between each entropy production component σand |X| is a single quadratic curve, (Fig. 2b). If the entropy production can be described by the Fourier law as due to heat conduction only, it can be expressed as . Hence, we obtain σ = 1.82 −  (we respectively used the minimum and maximum values of temperature for calculation because the temperature oscillates). This relation is consistent with the fitting curve in Fig. 2b. The results indicating that the single curves describe these relationships are justified by the fact that the thermodynamic flux is defined on the basis of a linear function of the thermodynamic force. Thus, for the entropy production expressed as a function of X, the MEPP is inapplicable for specifying the thermodynamic states and for predicting the transition point from an ASF to an HTW.

Figure 2

(a) Relationship between the local thermodynamic flux J and the local thermodynamic force X. (b) Relationship between the entropy production σ and the local thermodynamic force X for an ASF and HTW. The solid line is fit using the quadratic curve, . It is in agreement with the theoretical curve, .

We must therefore query how the MEPP successfully predicted the transition points of the nonequilibrium states in the experiments reported in refs[11,16,17]. In those experiments, instead of a local thermodynamic force, a driving force adjusted by the energy supplied from the surroundings was used as a thermodynamic force; for example, the degree of supersaturation, the degree of supercooling, and a pressure gradient were used. The driving force, which is constant in time, serves to prevent the internal system state from achieving thermal equilibrium. As the driving force increases, the system state changes. Thus, the driving force serves as a measure of the nonequilibrium degree. For our system, we next analyze the system’s behavior using the thermodynamic variables as functions of the driving force, which is determined by the temperature gradient between the inner and outer walls, i.e., F = (1/Tc − 1/Th)/(Ro − Ri). Note that F has the same dimension as a thermodynamic force owing to heat conduction.

Figure 3 shows the relationships between the driving force and each component of the local thermodynamic flux. The flow pattern changes from an ASF to an HTW with a jump in |J| at a driving force of F = 70–80 × 10−6 m−1K−1, corresponding to ΔT ≈ 7–8 K. Quantitative analysis reveals that for an ASF, the first derivative of |J| remains constant at the value of (3.07 ± 0.00(4)) × 106 Jm−1s−1K−1, whereas for an HTW, the value is (74.0 ± 27.3) × 106 Jm−1s−1K−1. In each flow pattern, |J| can be expressed by two straight lines with significantly different slopes as functions of F. For J, regardless of whether the flow pattern is an ASF or HTW, this term increases monotonically with respect to F. The values can be expressed using the theoretical relation J = λTcTh(Ro − Ri)F/Rln(Ro/Ri), where R(=30mm) is the measurement point. Finally, J slightly changes its derivatives, although it is not as well determined as J. From a quantitative analysis conducted, we found that for J, the slope of its derivative of an ASF region decreases by 41.8% in comparison with that of an HTW region. Interestingly, J and J are in a direction perpendicular to F and they exhibit a nontrivial change with F. HTW occurs owing to the symmetry-breaking of heat conduction. The thermodynamic fluxes perpendicular to F may provide a clue to a signature of the change in the flow patterns. According to the MEPP, a jump in the thermodynamic flux is the signature of a first order transition[33]. Therefore, we can distinguish one nonequilibrium state from another in terms of the thermodynamic fluxes with respect to the driving force. As J can be expressed by the linear function of F, the absolute value of entropy production calculated from Eq. (1) can be fitted by the following quadratic function:where L and θ are phenomenological coefficients[6,11,33]. The entropy production in each state can be expressed by different quadratic curves, i.e., for an ASF and HTW. We can also fit the experimental values of J using calculated from the relation .

Figure 3

Effect of the driving force F on each component of the local thermodynamic flux J for an ASF and HTW: (a–c) J, J, and J, respectively. (a) The solid line represents J = λTcTh(Ro − Ri)F/rln(Ro/Ri). (b) The solid and broken lines are two regression lines calculated from the relation σ = JX (See the text for more details). Inset: Magnified view.

The relationship between σ and F is shown in Fig. 4a. In the low-thermodynamic-force region, the ASF entropy production curve lies above that of the HTW, whereas the opposite is true for the high-thermodynamic-force region. The intersection point of the two curves is at F = 76.2 × 10−6 m−1K−1, which corresponds to ΔT = 7.59K when the fitting range for the HTW entropy production is 8–11 K (ΔT = 8.15 K for the 8-14-K range). We can interpret the system behavior according to the MEPP. That is, the ASF appeared below the point of intersection because the entropy production for the ASF exceeded that for the HTW. Above the point of intersection, the entropy production for the HTW became greater. The system changed to maximize the entropy production induced by J with respect to F. Thus, the transition point is the intersection point of the two curves for the ASF and HTW. Note that all absolute values of σ described here are the same values as σ described in Fig. 2(b). The only difference between the two entropy productions is the variable expressed by the function of X or F. σ with respect to F can be described as the two different curves in each flow pattern, whereas σ with respect to X falls on a single curve. The function σ(F) permits the distinction between the two flow patterns. θ in Eq. (2) is an important factor because, in the case of a nonzero value, state transition occurs. The physical meaning of θ was interpreted as a correction term for conversion to the local driving force[6]. Our results shown in Fig. 3(a) rule out this interpretation. Thus, θ may describe the interference of the two irreversible processes or the degree of symmetry-breaking because for σ, θ = 0. Full interpretation of the physical meaning needs to be investigated further.

Figure 4

(a) Effect of the driving force F on the y-component of the entropy production σ for an ASF and HTW. The solid and broken curves represent the regression curves for an ASF and HTW using Eq. (2): y = 5.18 × 10−8x2, and y = 3.00 × 10−5(x − 73.0)2, respectively. The intersection point of the two curves for the ASF and HTW correspond to the transition point. Inset: Magnified view. (b) Time course of the temperature T near the prediction points. ΔT is the temperature difference between the inner and outer walls.

To verify the transition point determined by the MEPP, we computed the time course of the temperature at the measurement position by varying ΔT in intervals of 0.1 K (Fig. 4b). At ΔT = 7.4 K, the temperature remained constant during the measurement time span. At ΔT = 7.5 K, however, the temperature started to increase at a time of approximately 260 s. Eventually, temperature oscillation began. This result indicates that HTW behavior appeared at ΔT = 7.5 K, and this transition point is extremely close to the value calculated from the MEPP. Although the calculation using the MEPP was performed in 1-K intervals, the precision of the predicted value is extremely high (note that we obtain complete agreement on the transition point if we set the intervals to 0.5 K). The MEPP correctly predict the transition point only in the time span in which the MEPP analysis is conducted. In fact, for ΔT = 7.4 K, the temperature oscillated at a time of approximately 400 s, and the HTW appeared. When the analysis period was changed from 200–300 s to 300–400 s, we found that the prediction point calculated by the MEPP decreased by approximately 0.1 K. These results indicate that, although the nonequilibrium system does not reach a steady state, the MEPP can specify the nonequilibrium state and predict the transition point.

This understanding motivates the question of how J is capable of uniquely specifying the nonequilibrium states. The decisive component of the thermodynamic variables that specifies the nonequilibrium states may be the component that breaks the spatial symmetry with respect to the driving force applied to the entire system. J at the measurement point is perpendicular to the driving force. In context, using J, we must be able to predict the transition point because J is perpendicular to the driving force. Therefore, using the same approach as that for the calculation of σ, we analyzed the transition points obtained from the σ component at different measurement points. For simplicity, we fit σ using a quadratic function of F, even though J may behave as a high-order dependency of F. The result is shown in Fig. 5.σ exhibits the distinctive features exhibited by σ: first, the system changes to maximize the entropy production; and second, the point of state transitions coincides with the intersection between the entropy production curves. These transition points are congruent with those obtained from the σ component (Fig. 5b).

Figure 5

(a) Entropy production for the z components, σ, as functions of the driving force F at various measurement points. The solid and broken lines represent the fitting curves of the entropy production for an ASF and HTW, respectively, obtained using Eq. (2). The intersection point of the two curves for the ASF and HTW correspond to the transition point. (b) Comparison of the transition-point prediction based on analysis of σ and σ at various measurement points.

We investigated the relationship between the driving force F and the total entropy production (Fig. 6). The total entropy production falls on a single curve, , which represents the theoretical value of the x component of entropy production expressed as a function of F. We unsuccessfully attempted to predict the transition point using the total entropy production as a function of F because the total entropy production increases monotonically with F. We found that information related to the system was lost because the σ component, which constitutes 94.6–99.9% of the total entropy production, covers the slight change in the entropy production due to the transition. This finding indicates that the largest component of the entropy production as well as the total entropy production cannot distinguish one nonequilibrium state from another. The results indicate that the current formulation of the MEPP requires revision.

Figure 6

Relationship between the total entropy production and the driving force F. The solid curve represents σ = λTT[(R − R)F/Rln(R/R)]2.

From the above results, we deduced that a specific component of the thermodynamic flux depending on the driving force distinguishes one nonequilibrium state from another. Hence, we may extend the MEPP using the relation between the thermodynamic flux and the driving force. This relation specifies the nonequilibrium states as thermodynamic phases similar to phases in equilibrium. Therefore, we can derive an equation related to the phase boundaries between the nonequilibrium states. Using the corrected thermodynamic force I = (F − θ), instead of the local thermodynamic force[6,11,33] and its conjugate variable J(F), the maximum entropy production principle can be rewritten aswhere μ is a Lagrangian multiplier. We obtain the following explicit expression for the thermodynamic flux:

Coefficient r corresponds to the degree of a homogeneous function[5]. If two different driving forces, F1, F2, are applied to a system, the entropy production can be expressed in the form σ = σ(F1, F2). Then, the differential form of the entropy production can be given as

This equation holds for different nonequilibrium phases. When two nonequilibrium phases are at a transition point, the entropy productions of the two phases are equal. Therefore, from Eqs (4) and (5), we obtainwhere ΔJ is the difference in the thermodynamic fluxes for the transition from one phase to another. This expression may correspond to the Clausius–Clapeyron equation describing the slope of the phase boundaries in equilibrium.

In this paper, the relationship between the thermodynamic flux and the various nonequilibrium states in a system incorporating two irreversible processes were examined. A nonequilibrium state changes to another state to maximize the entropy production induced by the thermodynamic flux in a direction perpendicular to the driving force. Hence, a selection rule for predicting the transition point between different states is needed to evaluate the symmetry-breaking component of the thermodynamic flux with respect to the driving force. Moreover, our analysis revealed that the total entropy production and the largest entropy production cover the slight change in the entropy production due to the transition. Thus, the quantities involved in this principle are likely to be system-dependent.

One wonders if the MEPP-based analysis can apply to nonequilibrium systems involving several branches of solutions, for example, bistability in flow patterns[3,34]. A slight difference in initial perturbation leads to a completely different patterns. In this manuscript, we used the steady value of driving force applied to the entire system as a kind of state variable describing entropy production and thermodynamic flux. This approach cannot be used for nonequilibrium systems exhibiting several branches of solutions because the same value of driving force may correspond to different entropy productions exhibited by different branches. We may use driving force reflecting initial conditions as state variables to predict the system behavior; the difference in the initial perturbation corresponds to the difference in the driving force. Therefore, when we use an initial value of driving force or its average value as nonlocal thermodynamic force, different values of the nonlocal thermodynamic force may specify different branches of solutions. If all branches of solutions show the same entropy production, the approach predicts that coexisting states will appear. If there is a slight difference in entropy production among several branches of solutions, the branch exhibiting the maximum entropy production will be selected. Consequently, applying our modified MEEP to nonequilibrium systems involving several branches of solutions is a significant issue to be addressed. Finally, we proposed an equation related to the phase boundaries between the nonequilibrium states on the basis of the MEPP involving entropy production as a function of the driving force. Experimental and numerical validation of this equation will reveal the range of MEPP applicability.

Numerical calculation methods

We investigate the thermal convection in an annular pool filled with silicon melt (height: 3 mm), for which there is a fixed temperature difference between the inner (Ri = 15 mm) and outer walls (Ro = 50 mm). The silicon melt is a noncompressible Newtonian fluid. We used the continuity equation as the governing equation along with the Navier–Stokes equation and energy equation:where is the fluid velocity vector, t is the time, ρ is the density of the silicone melt, is the pressure, ν is the kinetic viscosity of the silicone melt, T is the temperature, and α is the thermal diffusivity of the silicon melt. According to order analysis, the strength of natural convection relative to Marangoni convection can be determined by the parameter Ra1/2/Ma2/3(<1, Ra: Rayleighnumber, Ma: Marangoninumber). In this system, this value is small; therefore we ignored the effect of the gravity term in the Navier–Stokes equation. The boundary conditions at the free surface and at the container bottom are expressed by the following equations:

For the boundary condition at the free surface, we considered the generation of convection due to the thermal Marangoni effect. We also assumed that the shape of the free surface would not change. Solid walls were used for the inner wall, outer wall, and bottom surface of the container, and we used a no-slip condition for the velocity. We fixed the temperature of the inner wall to be constant and set the outer wall to a temperature higher than that of the inner wall. Specifically, we set the temperature of the inner wall to T = T = 1683 K (T: melting point of silicon), and we varied the temperature of the outer wall in the range of T(=ΔT + Tc) = 1684 K to 1697 K. We used adiabatic conditions for the upper surface and bottom surface. For the initial conditions, we set the velocity of the fluid within the container to 0, and the temperature was homogeneous at a value of 1683 K.

We set the number of grid points to be 81 in the radial direction, 180 in the circumferential direction, and 21 in the vertical direction. These numbers for the grid are based on the conditions used by Li et al.[35], and the resolution was sufficient to handle hydrothermal waves for a fluid with a low Prandtl number. We discretized the equations with the finite volume method and performed the numerical calculations using OpenFOAM, which uses the PISO algorithm. The details of the calculation method have been described in previous research[36]. The simulation code has already been verified, and we verified that there were no problems. The physical values used in the calculations were the same as the values used previously, as presented in Tables 1 and 2.

Table 1

Physical properties.

Physical properties	Value
Thermal conductivity (λ)	64 Wm−1K−1
Viscosity (μ)	7.0 × 10−4 kgm−1s−1
Density (ρ)	2530 kgm−3
Thermal expansion coefficient (β)	1.5 × 10−4 K−1
Surface tension coefficient (γT)	−7.0 × 10−5 Nm−1K−1
Heat capacity (Cp)	1000 Jkg−1K−1
Melting temperature (Tm)	1683 K

Table 2

Numerical parameters of the simulation.

Parameter	Value
Grashof number (Gr = gρ2βΔTd3/μ2)	1.19 × 104
Prandtl number (Pr = μCp/λ)	1.09 × 10−2
Marangoni number (Ma = −γTΔT(Ro − Ri)/μα, α = λ/ρCp)	5.54 × 102 − 2.08 × 103
Rayleigh number (Ra = Gr × Pr)	1.19 × 102