Analyzing Hydrothermal Waves Transitions through Rotational Field Application Based on Entropy Production

1. Introduction

Understanding the phenomena of transitions in the spatiotemporal structure of spontaneous convection under an applied external field is crucial in fundamental physics and industries. In fundamental physics, it relates to key areas such as nonlinear dynamics and nonequilibrium thermodynamics related to fluid instabilities [1,2,3,4]. In industries, these transitions are crucial for manufacturing process using transport phenomena [5,6,7,8,9]. Rayleigh–Bénard convection—a classic example of thermal convection—has been extensively studied both theoretically and experimentally for many years [10,11,12,13,14,15,16,17]. Chandrasekhar stated that the symmetry breaking of Bénard convection is sustained by the competition between the viscosity-induced dissipation of kinetic energy and the internal energy released by buoyancy [18]. He obtained stable conditions for convection generation. Hohenberg and Swift analyzed the complex behavior of Bénard convection under temporally modulated temperature differences using the Lorenz model [19]. They obtained a bifurcation diagram that illustrates the transition process not only from conductive heat transfer to Bénard convection but also from a hexagonal pattern to a roll pattern.

Recently, attempts have been made to theoretically analyze the evolution of nonequilibrium states in systems where various irreversible processes are coupled, using the variational principle of maximum entropy production (MEPP) [20,21,22,23,24]. The MEPP shows that nonequilibrium states evolve into more complex states that generate higher entropy. Here, "complex" refers to the interference between two or more irreversible processes [25]. When two fluids with different viscosities are displaced in a Hele–Shaw cell, the interface becomes destabilized, transitioning from a stable diffusion interface to a finger-like interface due to the Saffman–Taylor instability as the viscosity ratio increases [26]. With a further increase in the viscosity ratio, the flow becomes more complex owing to the splitting of the viscous fingers [27,28]. The flow at the growing interface evolves such that entropy production is maximized [29]. The transition of Bénard convection from a hexagonal to a rolled pattern owing to temperature modulation transitions to a pattern with higher entropy production [30]. The bistability of the hexagonal and rolled patterns arises because, while the rolled pattern exhibits higher heat flux, the hexagonal pattern still produces higher entropy.

In addition to predicting the state transitions, the variational principle can be used to derive the governing equations for complex transport phenomena. Doi et al. derived the time evolution equations for the contact line and droplet shape in the dynamics of an evaporating droplet, driven by the Marangoni effect [31]. They also derived time-evolution equations for deposition patterns, such as coffee rings and mountain shapes, during the drying process of the liquid film [32]. Thus, the MEPP has significantly advanced our understanding of various nonequilibrium phenomena. However, it remains unknown whether MEPP can be applied to complex transport phenomena involving multiple external forces. In this study, we quantitatively analyzed the state transitions of thermal convection in a system under a rotational field using the MEPP.

2. Equations of Motion and Numerical Details

Figure 1 shows the schematic of an annular pool with an applied temperature difference. A shallow annular pool with an inner wall of 15 mm $\left(R_{\mathrm{i}}\right)$, an outer wall of 50 mm $\left(R_{\mathrm{o}}\right)$, and a depth of 3 mm$\left(b\right)$ was filled with a silicon melt. The inner wall is fixed at the melting point of silicon, $T_c(=1683 \mathrm{K})$, and the outer wall is set at a higher temperature than the inner wall, $T_h(=T_c+\mathsf{\Delta }T)$. The upper liquid surface is a free interface, where Marangoni convection is used. No-slip conditions were imposed on the fluid velocities on all solid walls (inner and outer walls and bottom). The numerical analysis was performed in OpenFOAM using the PISO method, with discretization achieved through the finite volume method. The grid numbers are 81, 180, and 21 in the radial, circumferential, and vertical directions, respectively. This grid configuration, based on the conditions established by Li et al.,[33] sufficiently captures HTW in low-Pr fluids.

Assuming that the silicon melt behaves as an incompressible Newtonian fluid and the Boussinesq approximation is applicable, the continuity, Navier–Stokes, and energy equations can be expressed as

$${\displaystyle \frac{\partial v_i}{\partial x_i}}=0,$$

(1)

$${\displaystyle \frac{\partial v_i}{\partial t}}+v_i·{\displaystyle \frac{\partial v}{\partial x_i}}=-{\displaystyle \frac{1}{\rho }}{\displaystyle \frac{\partial p}{\partial x_i}}+\nu {\nabla }^2{v}_i,$$

(2)

$${\displaystyle \frac{\partial T}{\partial t}}+v_i·{\displaystyle \frac{\partial T}{\partial x_i}}=\alpha {\nabla }^{2}T,$$

(3)

where $v_i, \rho ,\nu ,$and$\alpha$are the velocity, density, kinematic viscosity, and thermal diffusivity of the fluid, respectively. Additionally, $p$ is pressure, $T$ is temperature, and $t$ is time. We assume that if a subscript appears more than once in the same term in the equation, summation over the repeated subscripts is implied. In this thermally driven system, spontaneous convection, resulting from the density difference, and Marangoni convection, resulting from the surface tension difference, coexist. However, the former can be neglected because the density changes owing to temperature differences are small. Therefore, it is unnecessary to consider the buoyancy term in Eq (2) [34]. Initially, there was no convection in the annular pool, and the fluid temperature was set uniformly at 1683 K, which is the melting point of silicon. Convection is induced by the Marangoni effect owing to the presence of a thermal gradient on the free surface, according to the following equation:

$$\eta {\displaystyle \frac{\partial v_x}{\partial z}}={\gamma }_{\mathrm{T}}{\displaystyle \frac{\partial T}{\partial x}},\eta {\displaystyle \frac{\partial v_y}{\partial z}}={\gamma }_{\mathrm{T}}{\displaystyle \frac{\partial T}{\partial y}},$$

(4)

where $\eta$ and ${\gamma }_{\mathrm{T}}$ are the viscosity and the surface tension coefficient of the fluid, respectively. The fluid properties and computational parameters used in the numerical calculations are listed in Table 1.

Increasing the temperature difference between the inner and outer walls results in the transition of the flow state in the system from an axially symmetric flow (ASF) to a rotating oscillatory flow called hydrothermal waves (HTW) [35]. The former is a radially symmetric flow, while the latter causes convection in the circumferential direction. The transition point of the flow state in the system is at approximately $\mathsf{\Delta }T=7 \mathrm{K}$. Figure 1b shows a typical ASF and HTW displayed in terms of surface temperature fluctuation $\delta T=T-1/2\pi {\int }_0^{2\pi }T\mathrm{d}\theta$. For ASF, the distribution is almost uniform, showing no significant temperature difference. In contrast, for HTW, a radial region with alternating hot and cold regions appears and rotates in the circumferential direction.

In a fluid undergoing an irreversible process, the rate of entropy production, $\sigma$, can expressed via the entropy produced per unit volume (s) and the conserved quantity (${\alpha }_k$) using the following equation [36].

$$\sigma ={\displaystyle \frac{\partial s}{\partial t}}={\displaystyle \frac{\partial s}{\partial {\alpha }_k}}{\displaystyle \frac{\mathrm{d}{\alpha }_k}{\mathrm{d}t}},$$

(5)

$$\eta {\alpha }_k=\left\{\begin{array}{c}\epsilon :\mathrm{energy}\\ \rho v :\mathrm{momentum} \\ {\rho }_m :\mathrm{mass}\end{array}\right\}$$

Subscript k represents the conserved quantity. Energy $\epsilon$ is transported by heat conduction, momentum $\rho v$ is transported by fluid movement, and matter ${\rho }_m$ is transported by mass transfer. Entropy is produced in the fluid because of the irreversible transport of conserved quantities through the driving force. Thus, the entropy production can be expressed as:

$$X_k={\displaystyle \frac{\partial s}{\partial {\alpha }_k}},J_k={\displaystyle \frac{\mathrm{d}{\alpha }_k}{\mathrm{d}t}}$$

(6)

$$\sigma =X_k{J}_k.$$

(7)

Here, $X_k$ represents the local driving force, and $J_k$ represents the transport of conserved quantities in response to the driving force. In thermal convection, two irreversible processes interact: the transport of energy via heat conduction and the transport of momentum via fluid movement. Because the entropy production of the latter is significantly smaller than that of the former, only the entropy production by the heat flux is considered here [34]. Therefore, subscript k is omitted. Experiments have demonstrated that even with complex flow patterns, such as those of a Bénard cell, the heat flux is linear with respect to the driving force [37,38]. This leads to the following equation:

$$X_i={\displaystyle \frac{\partial }{\partial x_i}}\left({\displaystyle \frac{1}{T}}\right),\hspace{1em}\hspace{1em}{J}_i=\lambda {\displaystyle \frac{\partial T}{\partial x_i}}$$

(8)

$$\sigma =X_i{J}_i$$

(9)

where $\lambda$ is the thermal conductivity of the fluid. Entropy production is a scalar quantity that represents the amount of heat produced per unit volume due to the thermal gradient in each direction. However, because certain components are highly related to the symmetry with respect to the driving force that keeps the system nonequilibrium, we will denote entropy production separately for each gradient in this study.

$${\sigma }_x=X_x{J}_x,{\sigma }_y=X_y{J}_y,{\sigma }_z=X_z{J}_z$$

(10)

As the system is far from equilibrium, both $X_i$ and $J_i$ are time-dependent and spatially distributed. To characterize the nonequilibrium state of the system through space- and time-dependent thermodynamic quantities, local thermodynamic quantities at point A = (30 mm, 0 mm, 3 mm) were calculated by time averaging.

$$\left|X_i\right|=1/\mathsf{\Delta }t{\int }_{t_o}^{t_1}\sqrt{{\left({\displaystyle \frac{\partial }{\partial x_i}}\left({\displaystyle \frac{1}{T}}\right)\right)}^2}\mathrm{d}t,\left|J_i\right|=1/\mathsf{\Delta }t{\int }_{t_o}^{t_1}\sqrt{{\left(\lambda {\displaystyle \frac{\partial T}{\partial x_i}}\right)}^2}\mathrm{d}t,$$

(11)

where $t_o$ = 200 s and $t_1$ = 300 s. Additionally, $\mathsf{\Delta }t$ represents the time interval $t_1-t_o$. The entropy production was calculated as the product of $\left|X_i\right|$ and $\left|J_i\right|$, as shown in Equation (9). In addition to the thermal gradient, a counterclockwise rotational field with $\omega =0, 0.5, 1, 3 {\mathrm{m}\mathrm{i}\mathrm{n}}^{-1}$ was applied to the annular pool to assess changes in the flow state and to measure heat flux and entropy production.

3. Results and Discussion

Figure 2 shows changes in temperature at point A on the fluid surface for different rotational velocities at $\mathsf{\Delta }T=14 \mathrm{K}$. In the absence of a rotating field, HTW occurs with a fluctuation range of approximately 0.7 K, with a period of approximately 7 s. For $\omega =0.5 {\mathrm{m}\mathrm{i}\mathrm{n}}^{-1}$, the HTW is a combination of long-term oscillations, with a period of approximately 90 s, and short-term oscillations, with a period of approximately 2 s. The temperature fluctuation of the long-term oscillation was comparable to that without rotation, whereas that of the short-term oscillation was approximately one-third that of the long-term oscillation. At $\omega =1.0 {\mathrm{m}\mathrm{i}\mathrm{n}}^{-1}$, the long-term oscillation disappeared, and the temperature fluctuation resulted in a constant-amplitude HTW with a period of approximately 1.5 s. For $\omega =3.0 {\mathrm{m}\mathrm{i}\mathrm{n}}^{-1}$, the temperature fluctuation disappeared, and the ASF appeared. Theoretical studies have shown that the addition of a rotating field stabilizes thermal convection [18].

Figure 3 shows the entropy production ${\sigma }_i$ in each direction calculated using Equations (10) and (11) as a function of the thermodynamic force $X_i$ at $\mathit{\omega }\mathbf{=}\mathbf{0.5}\mathbf{ }{\mathbf{m}\mathbf{i}\mathbf{n}}^{\mathbf{-}\mathbf{1}}$. As $X_i$ increases, ${\sigma }_i$ increases parabolically. The largest entropy production, ${\sigma }_x$, occurs in the x-direction due to heat flux, as it aligns with the thermal gradient on the inner and outer walls, which serves as the driving force. This is followed by ${\sigma }_y$ and ${\sigma }_z$. If the flow effect is ignored and entropy production is considered to solely result from heat conduction according to Fourier’s law, it can theoretically be expressed as ${\sigma }_i=\lambda T^2{X}_i^2$ [34]. When the theoretical values were plotted using the average temperature $T_{av}=\left(T_c+T_h\right)/2$, the acquired data closely matched the theoretical predictions. Even when the flow state changed from ASF to HTW, the entropy production data exhibited the same curve. This result accurately reflects the increase in entropy production owing to local driving forces but also indicates that changes in state cannot be distinguished by expressing entropy production solely as a function of thermodynamic forces.

Determining function to changes in nonequilibrium states is a critical issue. Previous experimental and computational examples have shown that using a function of the global driving force that maintains the nonequilibrium state of the system is effective [39,40]. The emergence of new flow states owing to symmetry breaking in response to the global driving force is often observed in various transport phenomena, such as Bénard convection, Marangoni convection, and the Saffman–Taylor instability [29,30,34,41]. A common feature of this phenomenon is the production of excess entropy in the direction orthogonal to the global driving force. This excess entropy is produced by thermodynamic flux in the direction of symmetry breaking. This thermodynamic flux is referred to as the emergent flux. When the emergent flux is considered as a function of the global driving force, the entropy production similarly becomes a function of the global driving force.

$$J_i=L\mathrm{’}\left(X_G-\mathsf{\Theta }\right)$$

(12)

$$\sigma =L{\left(X_G-\mathsf{\Theta }\right)}^2$$

(13)

where $L^{\prime }, L$, and $\mathsf{\Theta }$ are phenomenological coefficients. In particular, $\mathsf{\Theta }$ is significant as it is related to the free energy required to transition from one state to another and is also related to the degree of coherence of two or more irreversible processes [39,42]. Because the phenomenological coefficients are different for the different states of the system, the states of the system can be described by different lines. Using this relationship between the emergent flux and the global driving force, we can identify the state of the system. The entropy production for each nonequilibrium state can then be obtained using Equation (13). According to the MEPP, the thermodynamic flux changes such that the system transitions to a state with a higher entropy production. For the Bénard convection, the global driving force is the thermal gradient provided above and below the vessel[30,37]. As the driving force increased, the thermodynamic flux changed discontinuously with each transition from conduction to convective heat transfer and from hexagonal to roll patterns. Discontinuity represents the transition point of the nonequilibrium state, which is represented by the intersection of the entropy production curves. The MEPP, using the relationship between the emergent flux and global driving force, was used to quantitatively describe complex transitions in thermal convection under a rotating field. Additionally, it was used to predict the transition point.

In this thermally driven system, the thermal gradient provided on the inner and outer walls is the driving force that changes the state of the thermal convection inside. We express the thermal gradient in terms of the dimensions of the thermodynamic force, referring to it as the global driving force $X_G$.

$$X_G={\displaystyle \frac{1}{R_{\mathrm{o}}-R_{\mathrm{i}}}}\left({\displaystyle \frac{1}{T_c}}-{\displaystyle \frac{1}{T_h}}\right)$$

(14)

The heat fluxes in the x- and y-directions at each rotational speed are shown in Figure 4. The lower and upper abscissae represent the global driving force and the corresponding temperature difference $\mathsf{\Delta }T$, respectively. The numerical values of $\mathsf{\Delta }T$ are similar to those of $X_G$, except for an order of magnitude. $J_x$ increased linearly with respect to the global driving force regardless of the rotational speed. If the flow effects are disregarded and only heat conduction according to Fourier’s law is considered, the heat flux is theoretically expressed as $J_x=\lambda T_c{T}_h\left(R_{\mathrm{o}}-R_{\mathrm{i}}\right)X_G/\left\{r\ln \left(R_{\mathrm{o}}/R_{\mathrm{i}}\right)\right\}$[34]. The heat flux in the x-direction was almost consistent with the theoretical equation and was independent of the rotational speed. As the direction of the global thrust was radial, it agreed with the thermal gradient in the x-direction at the measurement point. Marangoni convection occurred in the x-direction at the free surface; however, the effect of the flow on the heat flux in the x-direction was negligible because the conduction heat transfer was more dominant. The heat flux in the y-direction, where the symmetry was broken by the global driving force, exhibited a nontrivial trend. In the absence of a rotational field, $J_y$ can be described by different straight lines for ASF and HTW using Equation (12). For ASF, $J_y$ passes through the origin and $\mathsf{\Theta }=0$, indicating a complete absence of interfering nature with other irreversible processes. For HTW, $\mathsf{\Theta }$ exhibits a finite value, indicating the degree of interference from other irreversible processes. In the presence of a rotational field, $J_y$ was on the same line as ASF, independent of the rotational speed; however, the value of the global driving force required to generate HTW increased significantly with increasing rotational speed. This indicates that the centrifugal force generated by the rotational field suppresses HTW generation. As the global driving force further increased, a discontinuous change in the heat flux was observed, and the slope of the straight line became even larger. This implies the existence of a different convection pattern.

The emergent heat flux results suggest the presence of a new HTW flow pattern. Figure 5 illustrates the surface temperature fluctuations for $\mathsf{\Delta }T=15 \mathrm{K}$ (top panels) and $\mathsf{\Delta }T=17 \mathrm{K}$ (bottom panels) at $\omega =1.0 {\mathrm{m}\mathrm{i}\mathrm{n}}^{-1}$. For $\mathsf{\Delta }T=17 \mathrm{K}$, time progresses from the left to the right panels at an interval of 0.6 s. Waves consisting of eight pairs of high- and low-temperature regions are generated from the inner wall, expanding radially while rotating counterclockwise. When the tip of the growing wave reaches half of the annular pool, the wave detaches from the inner wall and diffuses toward the outer wall. New waves reappear from the inner wall. This results in heat waves propagating in the circumferential direction, with waves being repeatedly generated and annihilated. A normal rotating oscillating flow with constant width and length that rotates in one direction is represented by HTW$\alpha$. On the other hand, a rotating oscillating flow that repeatedly generates and annihilates waves is represented as HTW$\beta$. In the presence of a rotating field, the heat flux region represented by the straight line with the steepest slope in Figure 4a corresponds to the region where HTW$\beta$ is observed.

The transitions from ASF to HTW$\alpha$ and from HTW$\alpha$ to HTW$\beta$ are discussed using MEPP. Figure 6 shows the emergent heat flux and the corresponding entropy production at $\omega =1.0 {\mathrm{m}\mathrm{i}\mathrm{n}}^{-1}$. To obtain the entropy production curve, we first determined the differential coefficient $\mathrm{d}{J}_y/\mathrm{d}{X}_G$ at each point of the heat flux. The heat flux was divided into three regions based on the differential coefficient values. In the regions divided by differential coefficients, the heat flux was fitted using Equation (12), providing three straight lines with different values of L’ and Θ. The entropy production in each region was fitted using Equation (13) with Θ in common, resulting in three entropy curves. The fitting parameters for each region are listed in Table 2.

The entropy production curves of ASF and HTW$\alpha$ intersect at $\mathsf{\Delta }{T}_c=12.4 \mathrm{K},$ while those of HTW$\alpha$ to HTW$\beta$intersect at $\mathsf{\Delta }{T}_c=16.3 \mathrm{K}$. This is nearly consistent with the point where the emergent heat flux changes discontinuously. Additionally, at $\omega =0.5 {\mathrm{m}\mathrm{i}\mathrm{n}}^{-1}$, the transition point was accurately predicted using MEPP. Three key observations can be observed from these results. First, MEPP effectively predicts the transition point of the state even in a nonequilibrium system with multiple external forces, including the addition of thermal gradients and rotational fields. Second, the relationship between the emergent heat flux and the global driving force leads to the discovery of new flow states. Third, the relationship between the emergent heat flux and global driving force allows the specification of nonequilibrium states to be specified as if they were a phase of an equilibrium system.

For example, in equilibrium, if the state of water can be categorized into solid, liquid, and vapor phases, each state can be represented by a distinct function, with the chemical potential $\mu$ expressed as a function of temperature (Figure 7a). The entropy per mole $S_m$ can be obtained from the slope as: ${\left(\partial \mu /\partial T\right)}_p=-S_m$. Because $S_m>0$, the chemical potential decreased as the temperature increased. The entropy increased in the order $S_{m.\mathrm{s}\mathrm{o}\mathrm{l}\mathrm{i}\mathrm{d}}<S_{m.\mathrm{l}\mathrm{i}\mathrm{q}\mathrm{u}\mathrm{i}\mathrm{d}}<S_{m.\mathrm{g}\mathrm{a}\mathrm{s}}$, leading to steeper gradients in the same order. For $T_{\mathrm{m}}<T$, extrapolation of $\mu$ for each phase shows that the solid phase is the most stable because it has the lowest μ. For $T_{\mathrm{m}}<T<T_{\mathrm{b}}$, the liquid phase is stable because it has the lowest the μ, while for $T>T_{\mathrm{b}}$, the gas phase is stable because it has the lowest μ. The phase transition point can be observed at the intersection of $\mu$ of each phase. These results are logical consequences of the principle of maximum entropy.

In nonequilibrium, if the emergent flux $J_{\perp }$ is expressed as a function of the global driving force $X_G$, each state can be described by a different function, allowing the nonequilibrium state to be defined as a phase, analogous to how chemical potential defines phases in equilibrium systems. As $X_G$ increases, $J_{\perp }$ increases, which can be accompanied by a discontinuity change at the transition point and an increase in slope. Therefore, the system evolves into a state with higher entropy production. Therefore, to achieve a higher entropy production rate, the system interferes with other irreversible processes and discontinuously changes into a state with higher flux. This leads to a state transition from ASF to HTW$\alpha$ and further to HTW$\beta$.

4. Conclusion

This study quantitatively analyzed the state transitions of the flow in a thermally convective system under an applied thermal gradient and rotational field using MEPP. Using a heat flux orthogonal to the applied thermal gradient, the three states of thermal convection, ASF, HTW$\alpha$, and HTW$\beta$, were described by distinct straight lines. The transition points between states are identified by the intersections of the entropy production curves for each state. This suggests that in nonequilibrium, state transitions are driven by changes in flux, as these changes lead to faster entropy production. The MEPP may have a certain universality, because it is found not only in thermal convection phenomena but also in various transport phenomena such as Marangoni convection and Saffman–Taylor instability. Because of the analogy with equilibrium thermodynamics in terms of phase description and stability, a theorem similar to the equilibrium system may be derived for the nonequilibrium system. Verification of the theorem will lead to further development of the MEPP.