Next Article in Journal
Response Time of Queueing Mechanisms
Previous Article in Journal
Soft Faint Continuity and Soft Faint Theta Omega Continuity between Soft Topological Spaces
Previous Article in Special Issue
Experimental Study on the Influence of Interfacial Energy Instability on the Flow Pattern Spatiotemporal Evolution of Thermal- Buoyant Capillary Convection
 
 
Font Type:
Arial Georgia Verdana
Font Size:
Aa Aa Aa
Line Spacing:
Column Width:
Background:
Article

Thermally Driven Convection Generated by Reaction Fronts in Viscous Fluids

Departamento de Ciencias, Sección Física, Pontificia Universidad Católica del Perú, Av. Universitaria 1801, Lima 15088, Peru
*
Author to whom correspondence should be addressed.
Symmetry 2024, 16(3), 269; https://doi.org/10.3390/sym16030269
Submission received: 30 January 2024 / Revised: 9 February 2024 / Accepted: 13 February 2024 / Published: 23 February 2024
(This article belongs to the Special Issue Fluid Flow and Heat Transfer, Symmetry and Asymmetry)

Abstract

:
Reaction fronts propagating in liquids separate reacted from unreacted fluid. These reactions may release heat, increasing the temperature of the propagating medium. As fronts propagate, they will induce density changes leading to convection. Exothermic fronts that propagate upward increase the temperature of the reacted fluid located underneath the front. For positive expansion coefficients, the warmer fluid will tend to rise due to buoyancy. In the opposite case, for fronts propagating downward with the warmer fluid on top, an unexpected thermally driven instability can also take place. In this work, we carry out a linear stability analysis introducing perturbations of fixed wavelength. We obtain a dispersion relation between the perturbation wave number and its growth rate. For either direction of propagation, we find that the front is stable for very short wavelengths, but is unstable for large enough wavelengths. We carry out a numerical solution of a cubic reaction–diffusion–advection equation coupled to Navier–Stokes hydrodynamics in a two-dimensional rectangular domain. We find transitions between the non-axisymmetric and axisymmetric fronts increasing with the width of the domain.

1. Introduction

Fluid motion may arise from buoyancy forces in liquids due to density changes. The density differences can be due to changes in the temperature of the liquid (such as in Rayleigh–Bénard convection) or changes in chemical composition. These mechanisms play an important role in diverse fields such as combustion [1,2,3], magnetohydrodynamics [4,5,6], and chemohydrodynamics [7,8]. The coupling of chemical reactions and hydrodynamics led to the observation of new instabilities and pattern formation in the presence of gravity. Experiments in the iodate–arsenous acid (IAA) reaction showed a transition to convection due to thermal and compositional changes across the front [9,10]. Experiments in the IAA reaction show flat reaction fronts propagating with constant speed in narrow tubes [11]. Increasing the diameter of the tube, the flat fronts became unstable evolving into curved fronts propagating with faster speeds. These fronts exhibit a transition from nonaxisymmetric to axisymmetric shapes as the diameter is further increased. Theoretical work accounted for the increase in speed and the front curvature in the limit of infinite thermal diffusivity [12]. Further experiments in Hele–Shaw cells exhibit finger splitting due to the combined thermal and compositional effects [13]. In the chlorite–tetrathionate (CT) reaction, thermal and compositional effects determine growth rates for finger formation in Hele–Shaw cells [14,15]. Further experiments for fronts in the CT reaction studied finger formation at different orientations [16]. Experiments in the CT reaction, carried out at T = 3 °C, were able to inhibit the thermal effects on the density changes [17]. Thermal and compositional effects on the density gradients determine the structures formed as the front propagates horizontally [18,19]. Fronts in the Fe(II)–Ni reaction propagating either upward or downward increase their speed due to thermal and compositional gradients [20]. Theoretical studies analyzed the hydrodynamic stability of thin fronts due to density gradients driven by thermal and compositional changes [21]. Fronts described by the Kuramoto–Sivashinski equation exhibit hydrodynamic instabilities generated by thermal effects [22]. Mukherjee and Paul analyzed fronts propagating in externally driven fluid flows [23,24]. They showed that thermal and compositional changes associated with the reaction increase the complexity of the front geometry.
In this manuscript, we consider reaction fronts separating reacted from unreacted fluids at different temperatures. Here, we focus on the thermal effects without including density changes due to chemical composition. The temperature difference is caused by heat released or absorbed in the reaction front. An exothermic reaction front propagating upward raises the temperature of the reacted fluid below the front. If the thermal expansion coefficient is positive, such as in the case of water at room temperature, the fluid underneath the front will have lower density. In this case, the less dense fluid will tend to rise, leading to a hydrodynamic instability. However, if the exothermic front propagates downward, with the less dense fluid above the front, an unexpected hydrodynamic instability may take place. Previous works studied this instability in the context of fluids confined in Hele–Shaw cells, modeling the hydrodynamics with Darcy’s Law [25,26]. They found the conditions of instability as a function of compositional and thermal gradients, and a Lewis number defined as the ratio between the thermal and molecular diffusivity. In this work, we consider effects of temperature on the stability of fronts traveling in viscous fluids. In this case, the stability of the front depends not only on a Lewis number and a Rayleigh number proportional to thermal density gradients. It also depends on a Schmidt number, defined as the ratio between the kinematic viscosity and the molecular diffusivity. In addition, we solve the nonlinear reaction–diffusion–advection equations for fronts propagating in a narrow rectangle resembling a two-dimensional tube. We found that near the instability, convection leads to fronts propagating with steady shape and constant velocity. The shape of these fronts can have different symmetries depending on the width of the domain and the initial conditions. There is a region of bistability where either type of front can take place.

2. Equations of Motion

We study the vertical propagation of reaction fronts in viscous fluids. The chemical front propagates along the vertical z direction, which points opposite to the gravity field. In the vertical direction there are no boundaries, so therefore the front can propagate indefinitely. We will set later the boundary conditions for the transverse direction. A reaction–diffusion–advection equation describes the time evolution of the front [27].
c t + V · c = D 2 c + f ( c ) ,
where c is the concentration of an autocatalytic reaction, V is the fluid velocity, D is the molecular diffusivity, and f ( c ) is the chemical reaction rate. In this work, we focus on the cubic autocatalytic reaction used to model the iodate-arsenous acid system [27,28]. For this reaction, we have f ( c ) = k c c 2 ( c 0 c ) , where c 0 is the concentration of the reacted fluid, and k c is a reaction rate constant. We use the Navier–Stokes equations in the Boussinesq approximation to describe viscous hydrodynamics. In this approximation, the effects of density variations are small except for in the large gravity term. Therefore, the continuity equation reduces to
· V = 0 ,
and the momentum equation becomes
V t + ( V · ) V = 1 ρ 0 P + ν 2 V ρ ρ 0 g z ^ .
Here, P is the pressure, ρ 0 is the mass density of the unreacted fluid initially at temperature T 0 , ρ is the fluid mass density, ν is the kinematic viscosity, g is the magnitude of the gravitational acceleration, and z ^ is a unit vector pointing along the z-axis. We consider only small amounts of heat involved in the reaction, leading to small density changes. As a result, we assume a linear relation between the fluid density ( ρ ) and temperature (T) given by
ρ = ρ 0 ( 1 α ( T T 0 ) ) ,
where α is the thermal expansion coefficient at constant pressure. The temperature follows an advection–diffusion equation that includes a heat source:
T t + V · T = D T 2 T + D T Q κ f ( c ) .
Here, D T represents the thermal diffusivity and κ represents the thermal conductivity. The parameter Q corresponds to the heat released by the reaction. We introduce the stream function ψ and the vorticity ω to solve this system for two-dimensional flow in the x-z plane. From Equation (2), the velocity components of the fluid can be expressed as V x = ψ / z and V z = ψ / x , and the vorticity is related to the stream function by ω = 2 ψ . Consequently, the Navier–Stokes equation [Equation (3)] can be written in terms of the vorticity and stream function
ω t = ψ x ω z ψ z ω x + ν 2 ω + g ρ 0 ρ x .
It is convenient to introduce units for time τ k c c 0 2 1 , and length D τ 1 / 2 , allowing us to define dimensionless coordinates x x / , t t / τ . We also introduce dimensionless variables ψ ψ 2 τ , ω ω τ , and c c / c 0 , T T / Δ T . Here, Δ T is the difference of temperatures between the reacted and unreacted fluid far away from the front. In addition, we define the dimensionless Schmidt number S c ν / D , Lewis number L e D T / D , and Rayleigh number
R a g α Δ T 3 ν D .
Therefore, the equations of motion in dimensionless form, after dropping the primes, are set as
c t = ψ x c z ψ z c x + 2 c + c 2 1 c ,
ω t = ψ x ω z ψ z ω x + S c 2 ω + R a S c T x ,
and
T t + V · T = L e 2 T + c 2 1 c .

2.1. Linear Stability Analysis

We carry out a linear stability analysis to determine the stability of the flat front solution, unbounded in the horizontal direction, by applying small perturbations to the concentration c = c ( 0 ) + c ( 1 ) . Here, c ( 0 ) corresponds to the concentration for the flat front solution without convection. Similarly, we introduce small perturbations to the temperature, the stream function, and the vorticity. We study perturbations of fixed wavenumber considering c ( 1 ) = c q ( z ) e σ t cos ( q x ) , T ( 1 ) = T q ( z ) e σ t cos ( q x ) , ψ ( 1 ) = ψ q ( z ) e σ t sin ( q x ) , and ω ( 1 ) = ω q ( z ) e σ t sin ( q x ) , where q and σ are the wavenumber and growth rate of the perturbations. We use a reference frame co-moving with the convectionless flat front traveling with constant speed v 0 in dimensionless units. Introducing this perturbation into Equation (8), keeping only the linear terms, and projecting over the corresponding cosine function, we obtain
σ c q = q ψ q d c ( 0 ) d z q 2 c q + d 2 c q d z 2 + c ( 0 ) 2 3 c ( 0 ) c q + v 0 d c q d z .
Similarly, by replacing the corresponding expression into the equation for the temperature Equation (10), we obtain
σ T q = q ψ q d T ( 0 ) d z L e q 2 T q + L e d 2 T q d z 2 + c ( 0 ) 2 3 c ( 0 ) c q + v 0 d T q d z ,
from Equation (9), we have
σ ω q = S c q 2 ω q + S c d 2 ω q d z 2 q R a S c T q + v 0 d ω q d z .
and from the relation between vorticity and stream function
ω q = d 2 ψ q d z 2 q 2 ψ q .
Given a wavenumber q, the linear system formed by Equations (11)–(14) leads to an eigenvalue equation for the growth rate σ with an infinite number of complex eigenvalues. The stationary flat front solution will be stable if all of the eigenvalues have a negative real part. Thus, to study the stability of the flat front solutions, we calculate the largest real part of σ for each wavenumber q.

2.2. Nonlinear Equations

To study the nonlinear front propagation, we consider a two-dimensional rectangular channel as a model of a vertical tube. The reaction starts at one end of the narrow rectangular domain and propagates along the length of the channel. We consider stress-free boundary conditions, so therefore the vorticity and the stream function are zero at the vertical walls. We set no flow boundary conditions at the walls for the temperature and the concentration ( c / x = 0 and T / x = 0 ). The temperature of the unreacted fluid far ahead of the front corresponds to T 0 = 0 , while far behind the front the temperature remains unchanged ( T / z = 0 ). In the channel boundaries, we impose stress-free boundary conditions, so therefore ψ = ω = 0. With these boundary conditions, we can compare the results of the linear stability analysis with the nonlinear calculations near the onset of convection.

3. Numerical Methods

We solve the two-dimensional Equations (8)–(10) using finite difference approximations for spatial derivatives over a rectangular mesh. The time evolution is carried out using the Euler method, considering the terms containing second derivatives implicitly, and the remaining terms explicitly. We program the algorithms using the FORTRAN programming language. We use a rectangular domain of length L z = 2600 and a variable width L x . The mesh size corresponds to Δ z = 0.5 in the vertical direction. We vary Δ x to test different widths, with the number of mesh points between 20 and 30 in the x direction. The time evolution is calculated using a first order finite difference approximation with temporal steps Δ t = 0.001 . Using this mesh size resulted in a flat front speed less than 0.3% different than the analytical result of 1/ 2 . We also tested our code for convective fronts with smaller spatial increments finding that our results do not change significantly. The front is initiated using a reacted concentration c = 1 in half the length and c = 0 in the other half. The front travels consuming the unreacted concentration c = 0 . As the front moves a certain distance, we shift all variables backwards, replacing the portion near the end of the domain with unreacted concentration ( c = 0 ) , zero stream function ( ψ = 0 ), and ambient temperature ( T = 0 ). This way the front can be kept away from the end of the channel. To obtain the growth rate for the linear Equations (11)–(13), we also use finite differences over a uniform mesh along the length L z . The equations are transformed into time evolution equations, which after a long time lead to the variables growing exponentially with the largest growth rate σ .

4. Results

4.1. Linear Stability Analysis

We study the stability of flat fronts calculating the dispersion relation σ versus q from Equations (11)–(13). The dispersion relation σ versus q corresponds to the growth rate with the largest real part for each wavenumber q. Since experiments take place in aqueous solutions, we can estimate the Schmidt number using the kinematic viscosity and the molecular diffusivity of water at the appropriate temperature [29,30]. Therefore, the Schmidt number can vary from 1246 at T = 4 °C to 138 at T = 50 °C. We set the Schmidt number to S c = 400 , which is close to experimental values for the iodate–arsenous acid systems [9]. This value can be considered constant since the change in temperature across the front is of the order of 0.1 °C. We use a Lewis number L e = 20 which is somewhat larger than typical values for the CT reaction ( L e = 10 ) [14,15,31], but smaller than values for the IAA reaction ( 50 L e 70 ) [9,32]. In the case of an exothermic reaction with a reaction front moving upwards, the warmer fluid is underneath the front. In this case, the Rayleigh number is positive if the coefficient of thermal expansion α is positive thus having the denser fluid on top. Figure 1 shows dispersion relations for positive Rayleigh numbers ( R a ) while keeping the Schmidt number constant ( S c = 400 ). In this figure, as we increase the wavenumber q away from q = 0 , the largest real part of the growth rate of the perturbations σ increases until it reaches a maximum positive value. It then decreases, becoming zero at q = q 0 , where it reaches critical stability. As we increase the Rayleigh number ( R a ), the value of the wavenumber at critical stability q 0 increases. The maximum value of the growth rate σ m a x also reaches higher values with increasing Rayleigh numbers. In this figure, we also show the location of the maximum growth rate on the dispersion relation curve ( q m a x , σ m a x ), where q m a x is the wavenumber corresponding to the maximum value of the growth rate σ m a x . This maximum on the dispersion relation curve shifts towards greater wavenumbers as the Rayleigh number increases. Therefore, the flat front becomes more unstable for larger values of R a , since the range of wavenumbers with positive growth rates increases. Figure 2 shows dispersion relations for different negative Rayleigh numbers ( R a ) having the same Schmidt number S c = 400 . As the magnitude of the Rayleigh number increases, the peak of the dispersion relation curves also shifts towards a greater wavenumber. The value of q 0 also increases, reflecting an increase in instability for larger absolute values of the Rayleigh number. We notice that the maximum values of σ are larger for negative values of R a ( Figure 1) than for positive values of R a ( Figure 2). The value of q 0 separates two regions of stability: perturbations with wavenumber q > q 0 will decrease with time, while perturbations for q < q 0 will increase exponentially. We show in Figure 3 the wavenumbers for critical stability q 0 as a function of Rayleigh numbers R a , while keeping the Schmidt number constant ( S c = 400 ). For either positive or negative Rayleigh numbers, the wavenumber q 0 increases with an increasing absolute value of R a . Therefore, the region of instability of the flat front solution becomes larger as the magnitude of the Rayleigh number increases. However, the increase in wavenumber q 0 is somewhat larger for negative Rayleigh numbers, compared to positive Rayleigh numbers of the same absolute value. We notice that for Rayleigh numbers close to zero, the value of q 0 is nonzero. Consequently, very small Rayleigh numbers (positive or negative) can destabilize the front. Previous studies show that this is the case for Darcy’s law hydrodynamics with positive Rayleigh numbers, whereas negative Rayleigh numbers can destabilize the front for values beyond a certain threshold, using smaller Lewis numbers [25].
The system of Equations (11)–(14) corresponds to eigenvalue equations for the growth rate σ and the corresponding eigenfunctions for the concentration c q , the temperature T q , the stream function ψ q , and the vorticity ω q . We continue to set the Schmidt number as S c = 400 . In Figure 4, we show the largest eigenvalues corresponding to the stream function ψ q at the critical value defined by q 0 . Here, we compare the results for R a = 2 and R a = 2 with the front centered at the position where c 0 = 0.5 . The eigenfunctions are normalized such that the maximum absolute value of ψ q is equal to 1. For R a = 2 , we notice a positive maximum value, followed by a minimum of negative value. This indicates convection forming above and below the front, with the strongest flow below the front. For R a = 2 , the situation is similar, however, the positive maximum below the front is smaller, and the negative minimum above the front is larger. Convection now will be stronger above the front. Changing the sign of the eigenfunctions only reverses the direction of convective flow, without modifying the strength of the convective flow. In Figure 5, we also compare the eigenfunctions at the wavenumber for maximum growth rate ( q m a x ) for R a = 2 and R a = 2 . In the case of the negative value ( R a = 2 ) , the stream function ψ q has a single maximum, located slightly above the front. In contrast, for R a = 2 there is a maximum and a minimum, indicating convection above and below the front. The eigenvalues for the stream function show that convection near the instability varies significantly depending on the sign of R a .
We study the dependence of the transition to convection for flat fronts on the Schmidt number. Figure 6 shows dispersion relations σ versus q with the same positive Rayleigh numbers, but different Schmidt numbers. Figure 6a,b corresponds to R a = 5 and R a = 10 , respectively. The curves for S c = 400 and S c = 4000 are close to each other, compared to the other curves. The values for the wavenumber at critical stability ( q 0 ) are almost the same in all curves, except for S c = 4 where it is somewhat larger. In each figure, we notice that larger Schmidt numbers lead to larger maximum values of the growth rate σ m a x . We also display the trajectory of q m a x as the Schmidt number varies. It shows that q m a x reaches a maximum value as a function of S c . This behavior is similar in both Figure 6a,b, but for a given value of q, the value of σ is larger for R a = 10 .
In Figure 7, we display dispersion relations for different Schmidt numbers, while keeping the same value of the negative Rayleigh number. Figure 7a,b corresponds to R a = 1 and R a = 2 , respectively. The curves are similar to each other and to the previous Figure 6 where the Rayleigh number is positive. The dispersion relations increase from zero at near q = 0 , reaching a maximum value at q m a x , then becoming negative after a critical wavenumber q 0 . Compared with positive Rayleigh numbers, the values of the maximum growth rate σ m a x are larger for negative Rayleigh numbers. This is significant when considering that the absolute values of R a are smaller in Figure 7. We also display the trajectory of the wavenumber for maximum growth rate, showing that it reaches a maximum value as the Schmidt number varies. While the wavenumber q m a x varies, the values for critical stability q 0 are very close in each figure. For R a = 1 , the value of q 0 is close to 0.2 , and for R a = 2 , its value is close to 0.25 . For each R a , the curves are very close to each other for large q, but the behavior is different near the origin. The dispersion relations yield higher values for σ if the Schmidt number increases. In contrast with Figure 6, the curves differ further as we increase S c from 400 to 4000. Although the curves have similar shape, their behavior as a function of S c is different for R a < 0 than for R a > 0 . This may indicate that approximating the hydrodynamic equations in the large S c limit may not lead to accurate results for R a < 0 .
As we increase the Schmidt number, the peaks on the dispersion relations shift slightly towards smaller values of the wavenumber q while the heights of the peaks increase considerably. However, the wavenumber corresponding to the critical stability value ( σ = 0 ) occurs near q = 0.2 , independently of the Schmidt number. For greater values of the wavenumber ( q > q 0 ), there is no quantitative difference in the dispersion relation obtained with different Schmidt numbers. When we increase the magnitude of the negative Rayleigh number ( R a = 2 ), the growth rate σ m a x of the dispersion relations reaches higher values as shown in Figure 7b, as compared to the previous case ( R a = 1 ). As we increase the Schmidt number, the peaks on the dispersion relation shift towards smaller values of the wavenumber q. However, the wavenumber corresponding to the critical stability value ( σ = 0 ) appears near q = 0.25 and remains at this value independently of the Schmidt number used to obtain the dispersion relation. The increment of the Schmidt number increases the slope of the dispersion relations near q = 0 . Thus, a high value of the Schmidt number contributes to destabilizing the flat front solutions when they are subject to small perturbations of long wavelength.
We study the effects of the Schmidt number on the linear stability analysis by calculating the maximum values of the growth rate σ m a x and the corresponding wavenumber ( q m a x ) . We analyze these values, keeping the Rayleigh number constant while varying the Schmidt number. In Figure 8, we show the values of q m a x as a function of the Schmidt number, whereby each curve corresponds to a different Rayleigh number. In all curves, as the Schmidt number increases from S c = 0.5 , the value of q m a x increases until it reaches a maximum. For positive Rayleigh numbers, the values of q m a x become almost constant for large values of S c . However, for negative Rayleigh numbers, q m a x continues to decrease significantly at larger S c . We notice that the curve for R a = 10 is above the curve for R a = 5 , indicating that perturbations with a shorter wavelength will grow faster if the Rayleigh number is higher. This is also the case for the curves with negative Rayleigh numbers, considering the absolute value of R a since the curve with R a = 2 is above the one with R a = 1 . We show in Figure 9 the dependence of σ m a x as a function of S c . The values of σ m a x for a given R a increase monotonically as a function of S c . However, for R a > 0 , the curves appear to reach a constant value for large S c . In contrast, the curves for R a < 0 increase more rapidly for larger S c . This result shows that using a large Schmidt number approximation (Stokes flow) may provide good results for positive R a . For negative R a , Figure 9 shows that the results continue to vary at large S c .

4.2. Nonlinear Front Propagation

We set a reaction front to propagate in a narrow rectangular domain resembling a two-dimensional tube. An exothermic front propagating upwards with a positive Rayleigh number has a less dense fluid underneath a denser fluid in regions away from the front. We find that for narrow widths, convection does not take place, and the flat front is a stable solution. As the width of the domain widens, the flat front loses stability. For stress-free boundary conditions, the results of the linear stability analysis can be used to obtain the critical width for the onset of convection. This width corresponds to L x = π / q 0 , where q 0 is the wavenumber for critical instability ( σ = 0 ). Increasing the width beyond this critical width, the front develops a constant shape, propagating with constant velocity. The speed of this front is higher than the flat front speed due to convection. Slightly above this critical width, convection evolves into two counterrotating convective rolls, one above the other. These rolls are steady in a reference frame traveling with the front (Figure 10). In this system, buoyancy drives the formation of convection. For each roll, the less dense fluid rises near a domain wall, with heavier fluid falling on the opposite wall. The resulting curved front is higher on one side of the channel. Although the lower roll is stronger, the front is affected by both, leading to a nonaxisymmetric shape. Here, we define the central axis as a vertical line equidistant from the walls. In Figure 11a, we show a nonaxisymmetric reaction front formed at a domain width slightly above the critical width. Figure 12 shows the speeds of fronts with steady shape as a function of domain width. We notice the sudden increase in speed at the transition from the flat to the nonaxisymmetric front. Increasing the width of the domain further, we find that the speed of the front increases until it reaches a maximum speed. As we continue to increase the width, the speed decreases, reaching a value where the nonaxisymmetric front loses stability. At this point, there is a transition to an axisymmetric front. The axisymmetric front is lifted higher on the axis, and lower at the walls, as shown in Figure 11b. We track the axisymmetric solutions using it as an initial condition, varying the width only slightly. By repeating this, we find axisymmetric solutions for different widths, up to a width of L x = 100 . Continuing towards smaller widths, the solution loses stability around L x = 48 , where the front makes a transition to the nonaxisymmetric branch. Figure 12 shows a region of bistability between the axisymmetric and nonaxisymmetric fronts. In this region, the speed of the axisymmetric fronts is higher, except for a small interval. This interval is near the intersection of the two branches in Figure 12. We point out that the speeds of the axisymmetric fronts are almost the same as the speeds of a nonaxisymmetric front evaluated at half the width. This indicates that the axisymmetric front corresponds to a nonaxisymmetric solution and its mirrored reflection, placed side by side, as suggested by Figure 12. However, the stability conditions for each branch are different.
We can characterize the strength of the convective flow with the maximum absolute value of the stream function ( | ψ | m a x ). Figure 13 displays this value as a function of width for the axisymmetric and nonaxisymmetric solutions. In each branch, | ψ | m a x increases monotonically, contrary to the speed (Figure 12) that reaches a maximum in each case. We also point out that the values of | ψ | m a x for the axisymmetric branch correspond to the values of the nonaxsimmetric branch evaluated at half the width.
We also calculate the speed of the front as a function of width for the negative Rayleigh number R a = 1.2 . In this case, we also find a nonaxisymmetric front just above the onset of convection, as shown in Figure 14. The speed increases monotonically, up to a width where the solution does not show a constant shape. In this branch, the solutions have higher speeds compared to the solution with the positive Rayleigh number of the same magnitude (Figure 12). In contrast, we did not find a region of bistability with an axisymmetric state. The axisymmetric solutions appear over a smaller width interval compared to Figure 12.

5. Discussion

We studied the hydrodynamic stability of reaction fronts in viscous fluids due to heat released or absorbed by the reaction. The fronts obey a reaction–diffusion–advection equation with cubic autocatalysis. The stability of the fronts propagating in viscous fluids depends on three dimensionless numbers: a Rayleigh number ( R a ), a Schmidt number ( S c ), and a Lewis number ( L e ). In this work, we kept the Lewis number equal to 20. The Rayleigh number can be either positive or negative depending on the type of heat transfer from the reaction, the direction of propagation, and the sign of the thermal expansion coefficient. The linear stability analysis shows that unbounded fronts will always be unstable, regardless of the Rayleigh number (positive or negative). For a given Rayleigh number, the dispersion relation σ versus q depends on the value of the Schmidt number. We observed that the maximum growth rate σ m a x increases with the increasing magnitude of R a . However, for higher Schmidt numbers, they become closer to each other for R a > 0 , while becoming further apart for R a < 0 . This shows the importance of the Schmidt number for front propagation. Numerical solutions of the full nonlinear equations showed a transition to convection for fronts propagating vertically in narrow rectangular domains. The transition takes place as the width of the domain increases. When the flat front becomes unstable, the system evolves into fronts of constant shape and speed. Near the transition, the fronts take a nonaxisymmetric shape, with one side of the front raised above the other. The speed increases as the width of the domain increases, until the front loses stability. For positive Rayleigh numbers, the front becomes a steady axisymmetric front. In this case, we found a region of bistability between the axisymmetric and nonaxisymmetric fronts. For negative Rayleigh numbers, the front becomes nonaxisymmetric near the flat front instability. Increasing the width leads to fronts with complex behavior. Further theoretical work may characterize these time-dependent states. Experiments in the IAA reaction [9,13], the CT reaction [15,17], and the Fe(II)–Ni reaction [20] determined the importance of thermally driven convection in propagating fronts. The present theory accounts for convection due to thermal effects. In the case of fronts propagating upwards in the IAA reaction, we find a critical wavelength ( λ = 2 π / q 0 ) equal to 0.36 cm, comparable to the value of 0.52 cm obtained using the eikonal relation [33]. These results validate the conclusion that compositional gradients drive the flat front instability, since experiments estimate a value of λ = 1.29 mm [9]. Another important conclusion from this work is that thermal effects drive convection in fronts propagating upward and downward. This conclusion is relevant to experiments in the CT reaction where removing thermal effects showed a reduced instability for propagating fronts in either direction [17]. A direct quantitative comparison with the experiments would require compositional effects, as well as heat losses in each particular geometry. Experiments can be designed where the exothermicity of the reaction is higher, such as in the case of combustion.

6. Conclusions

In summary, we found that upward and downward propagating reaction fronts in liquids become unstable due to heat released by the reaction. This leads to an unexpected instability where the less dense fluid is above the denser fluid. We showed that for large Schmidt numbers and positve Rayleigh numbers, the critical wavelengths approach a constant value. However, for negative Rayleigh numbers these values diverge. This theory is relevant to previous experiments with exothermic reactions. Further theoretical work should include the effects due to compositional gradients, as well as the particular geometry of the system.

Author Contributions

Conceptualization, P.M.V., R.G. and D.A.V.; methodology, P.M.V., R.G. and D.A.V.; software, P.M.V., R.G. and D.A.V.; validation, P.M.V., R.G. and D.A.V.; formal analysis, P.M.V., R.G. and D.A.V.; investigation, P.M.V., R.G. and D.A.V.; resources, P.M.V., R.G. and D.A.V.; data curation, P.M.V.; writing—original draft preparation, P.M.V. and D.A.V.; writing—review and editing, P.M.V., R.G. and D.A.V.; visualization, P.M.V. and D.A.V.; supervision, D.A.V. All authors have read and agreed to the published version of the manuscript.

Funding

R.G. acknowledges funding from FONDECYT through Grant No. 236-2015.

Data Availability Statement

The data that support the findings of this study are available from the corresponding author upon reasonable request.

Conflicts of Interest

The authors declare no conflicts of interest.

References

  1. Sivashinsky, G.I. Nonlinear analysis of hydrodynamic instability in laminar flames—I. Derivation of basic equations. Acta Astronaut. 1977, 4, 1177–1206. [Google Scholar] [CrossRef]
  2. Altantzis, C.; Frouzakis, C.E.; Tomboulides, A.G.; Matalon, M.; Boulouchos, K. Hydrodynamic and thermodiffusive instability effects on the evolution of laminar planar lean premixed hydrogen flames. J. Fluid Mech. 2012, 700, 329–361. [Google Scholar] [CrossRef]
  3. Matalon, M.; Metzener, P. The effect of thermal expansion on diffusion flame instabilities. J. Fluid Mech. 2010, 647, 453–472. [Google Scholar] [CrossRef]
  4. Fetecau, C.; Shah, N.A.; Vieru, D. General Solutions for Hydromagnetic Free Convection Flow over an Infinite Plate with Newtonian Heating, Mass Diffusion and Chemical Reaction. Commun. Theor. Phys. 2017, 68, 768–782. [Google Scholar] [CrossRef]
  5. Qayyum, S.; Hayat, T.; Alsaedi, A.; Ahmad, B. MHD nonlinear convective flow of thixotropic nanofluid with chemical reaction and Newtonian heat and mass conditions. Results Phys. 2017, 7, 2124–2133. [Google Scholar] [CrossRef]
  6. Ahmed, S.; Zueco, J.; López-González, L.M. Effects of chemical reaction, heat and mass transfer and viscous dissipation over a MHD flow in a vertical porous wall using perturbation method. Int. J. Heat Mass Transf. 2017, 104, 409–418. [Google Scholar] [CrossRef]
  7. De Wit, A. Chemo-Hydrodynamic Patterns and Instabilities. Annu. Rev. Fluid Mech. 2020, 52, 531–555. [Google Scholar] [CrossRef]
  8. De Wit, A.; Eckert, K.; Kalliadasis, S. Introduction to the Focus Issue: Chemo-Hydrodynamic Patterns and Instabilities. Chaos 2012, 22, 037101. [Google Scholar] [CrossRef] [PubMed]
  9. Edwards, B.F.; Wilder, J.W.; Showalter, K. Onset of convection for autocatalytic reaction fronts: Laterally unbounded system. Phys. Rev. A 1991, 43, 749–760. [Google Scholar] [CrossRef] [PubMed]
  10. Pojman, J.A.; Epstein, I.R.; McManus, T.J.; Showalter, K. Convective effects on chemical waves. 2. Simple convection in the iodate-arsenous acid system. J. Phys. Chem. 1991, 95, 1299–1306. [Google Scholar] [CrossRef]
  11. Masere, J.; Vasquez, D.A.; Edwards, B.F.; Wilder, J.W.; Showalter, K. Nonaxisymmetric and axisymmetric convection in propagating reaction-diffusion fronts. J. Phys. Chem. 1994, 98, 6505–6508. [Google Scholar] [CrossRef]
  12. Wilder, J.W.; Vasquez, D.A.; Edwards, B.F. Nonlinear front evolution of hydrodynamic chemical waves in vertical cylinders. Phys. Rev. E 1997, 56, 3016–3019. [Google Scholar] [CrossRef]
  13. Sebestíková, L.; D’Hernoncourt, J.; Hauser, M.J.B.; Müller, S.C.; De Wit, A. Flow-field development during finger splitting at an exothermic chemical reaction front. Phys. Rev. E 2007, 75, 026309. [Google Scholar] [CrossRef]
  14. Bánsági, T., Jr.; Horváth, D.; Tóth, A. Multicomponent convection in the chlorite–tetrathionate reaction. Chem. Phys. Lett. 2004, 384, 153–156. [Google Scholar] [CrossRef]
  15. Bánsági, T., Jr.; Horváth, D.; Tóth, A.; Yang, J.; Kalliadasis, S.; De Wit, A. Density fingering of an exothermic autocatalytic reaction. Phys. Rev. E 2003, 68, 055301. [Google Scholar] [CrossRef] [PubMed]
  16. Horváth, D.; Bánsági, T., Jr.; Tóth, A. Orientation-dependent density fingering in an acidity front. J. Chem. Phys. 2002, 117, 4399–4402. [Google Scholar] [CrossRef]
  17. Tóth, T.; Horváth, D.; Tóth, A. Thermal effects in the density fingering of the chlorite–tetrathionate reaction. Chem. Phys. Lett. 2007, 442, 289–292. [Google Scholar] [CrossRef]
  18. Rongy, L.; Schuszter, G.; Sinkó, Z.; Tóth, T.; Horváth, D.; Tóth, A.; De Wit, A. Influence of thermal effects on buoyancy-driven convection around autocatalytic chemical fronts propagating horizontally. Chaos 2009, 19, 023110. [Google Scholar] [CrossRef] [PubMed]
  19. Schuszter, G.; Pótári, G.; Horváth, D.; Tóth, A. Three-dimensional convection-driven fronts of the exothermic chlorite-tetrathionate reaction. Chaos 2015, 25, 064501. [Google Scholar] [CrossRef]
  20. Pojman, J.A.; Nagy, I.P.; Epstein, I.R. Convective effects on chemical waves. 3. Multicomponent convection in the iron (II)-nitric acid system. J. Phys. Chem. 1991, 95, 1306–1311. [Google Scholar] [CrossRef]
  21. Vasquez, D.A.; Edwards, B.F.; Wilder, J.W. Finite thermal diffusivity at onset of convection in autocatalytic systems: Discontinuous fluid density. Phys. Fluids 1995, 7, 2513–2515. [Google Scholar] [CrossRef]
  22. Guzman, R.; Vasquez, D.A. Front instabilities in the presence of convection due to thermal and compositional gradients. Chaos 2024, 34, 013123. [Google Scholar] [CrossRef] [PubMed]
  23. Mukherjee, S.; Paul, M.R. Velocity and geometry of propagating fronts in complex convective flow fields. Phys. Rev. E 2019, 99, 012213. [Google Scholar] [CrossRef] [PubMed]
  24. Mukherjee, S.; Paul, M.R. Propagating fronts in fluids with solutal feedback. Phys. Rev. E 2020, 101, 032214. [Google Scholar] [CrossRef] [PubMed]
  25. D’Hernoncourt, J.; De Wit, A.; Zebib, A. Double-diffusive instabilities of autocatalytic chemical fronts. J. Fluid Mech. 2007, 576, 445–456. [Google Scholar] [CrossRef]
  26. D’Hernoncourt, J.; Zebib, A.; De Wit, A. Reaction Driven Convection around a Stably Stratified Chemical Front. Phys. Rev. Lett. 2006, 96, 154501. [Google Scholar] [CrossRef]
  27. Leconte, M.; Martin, J.; Rakotomalala, N.; Salin, D.; Yortsos, Y.C. Mixing and reaction fronts in laminar flows. J. Chem. Phys. 2004, 120, 7314–7321. [Google Scholar] [CrossRef]
  28. Hanna, A.; Saul, A.; Showalter, K. Detailed studies of propagating fronts in the iodate oxidation of arsenous acid. J. Am. Chem. Soc. 1982, 104, 3838–3844. [Google Scholar] [CrossRef]
  29. Kestin, J.; Sokolov, M.; Wakeham, W.A. Viscosity of liquid water in the range -8 °C to 150 °C. J. Phys. Chem. Ref. Data 1978, 7, 941–948. [Google Scholar] [CrossRef]
  30. Holz, M.; Heil, S.R.; Saccob, A. Temperature-dependent self-diffusion coefficients of water and six selected molecular liquids for calibration in accurate H NMR PFG measurements. Phys. Chem. Chem. Phys. 2000, 2, 4740–4742. [Google Scholar] [CrossRef]
  31. Kalliadasis, S.; Yang, J.; De Wit, A. Fingering instabilities of exothermic reaction-diffusion fronts in porous media. Phys. Fluids 2004, 16, 1395–1409. [Google Scholar] [CrossRef]
  32. Gérard, T.; De Wit, A. Stability of exothermic autocatalytic fronts with regard to buoyancy-driven instabilities in presence of heat losses. Wave Motion 2011, 48, 814–823. [Google Scholar] [CrossRef]
  33. Wilder, J.W.; Edwards, B.F.; Vasquez, D.A. Finite Thermal Diffusivity at Onset of Convection in Autocatalytic Systems: Continuous Fluid Density. Phys. Rev. A 1992, 45, 2320–2327. [Google Scholar] [CrossRef]
Figure 1. Dispersion relation for different positive Rayleigh numbers ( R a ). The Schmidt number is set to S c = 400 . The orange curve displays the position of the maxima as the Rayleigh number varies.
Figure 1. Dispersion relation for different positive Rayleigh numbers ( R a ). The Schmidt number is set to S c = 400 . The orange curve displays the position of the maxima as the Rayleigh number varies.
Symmetry 16 00269 g001
Figure 2. Dispersion relation for different negative Rayleigh numbers ( R a ). The Schmidt number is set to S c = 400 . The orange curve displays the position of the maxima as the Rayleigh number varies.
Figure 2. Dispersion relation for different negative Rayleigh numbers ( R a ). The Schmidt number is set to S c = 400 . The orange curve displays the position of the maxima as the Rayleigh number varies.
Symmetry 16 00269 g002
Figure 3. The wavenumber q 0 for critical stability ( σ = 0 ) for positive and negative Rayleigh numbers ( R a ). The Schmidt number is set to S c = 400 .
Figure 3. The wavenumber q 0 for critical stability ( σ = 0 ) for positive and negative Rayleigh numbers ( R a ). The Schmidt number is set to S c = 400 .
Symmetry 16 00269 g003
Figure 4. Eigenvalues for the coefficients of the stream function ψ q for σ = 0 . Here, q 0 = 0.1785 for R a = 2 and q 0 = 0.2639 for R a = 2 . The eigenfunctions are scaled such that the maximum absolute value of each function is equal to 1. The Schmidt number is set to S c = 400 .
Figure 4. Eigenvalues for the coefficients of the stream function ψ q for σ = 0 . Here, q 0 = 0.1785 for R a = 2 and q 0 = 0.2639 for R a = 2 . The eigenfunctions are scaled such that the maximum absolute value of each function is equal to 1. The Schmidt number is set to S c = 400 .
Symmetry 16 00269 g004
Figure 5. Eigenvalues for the coefficients of the stream function ψ q at the maximum growth rate. We show two cases, one with R a = 2 ( q m a x = 0.0495 , σ = 0.27692 ) and another for R a = 2 ( q m a x = 0.0424 , σ = 2.183 ). The eigenvalue curves have been rescaled such that the peaks equal unity. The Schmidt number is set to S c = 400 .
Figure 5. Eigenvalues for the coefficients of the stream function ψ q at the maximum growth rate. We show two cases, one with R a = 2 ( q m a x = 0.0495 , σ = 0.27692 ) and another for R a = 2 ( q m a x = 0.0424 , σ = 2.183 ). The eigenvalue curves have been rescaled such that the peaks equal unity. The Schmidt number is set to S c = 400 .
Symmetry 16 00269 g005
Figure 6. The dispersion relation σ versus q for different values of the Schmidt number ( S c ). The orange curve displays the position of the maxima as the Schmidt number varies. (a) The Rayleigh number R a = 5 . (b) The Rayleigh number R a = 10 .
Figure 6. The dispersion relation σ versus q for different values of the Schmidt number ( S c ). The orange curve displays the position of the maxima as the Schmidt number varies. (a) The Rayleigh number R a = 5 . (b) The Rayleigh number R a = 10 .
Symmetry 16 00269 g006
Figure 7. The dispersion relation σ versus q for different values of the Schmidt number ( S c ). The orange curve displays the position of the maxima as the Schmidt number varies. (a) The Rayleigh number R a = 1 . (b) The Rayleigh number R a = 2 .
Figure 7. The dispersion relation σ versus q for different values of the Schmidt number ( S c ). The orange curve displays the position of the maxima as the Schmidt number varies. (a) The Rayleigh number R a = 1 . (b) The Rayleigh number R a = 2 .
Symmetry 16 00269 g007
Figure 8. Wavenumber q m a x corresponding to the maximum growth rate ( σ m a x ) on the dispersion relation as a function of the Schmidt number ( S c ) for different Rayleigh numbers ( R a ). In the case of negative Rayleigh numbers, q m a x continues to decrease for large values of S c , while for positive Rayleigh numbers, they approach a constant value.
Figure 8. Wavenumber q m a x corresponding to the maximum growth rate ( σ m a x ) on the dispersion relation as a function of the Schmidt number ( S c ) for different Rayleigh numbers ( R a ). In the case of negative Rayleigh numbers, q m a x continues to decrease for large values of S c , while for positive Rayleigh numbers, they approach a constant value.
Symmetry 16 00269 g008
Figure 9. The maximum value of the growth rate ( σ m a x ) on the dispersion relation as a function of the Schmidt number ( S c ) for different Rayleigh numbers ( R a ). In the case of negative Rayleigh numbers, σ m a x continues to increase for large values of S c . While for positive Rayleigh numbers, they approach a constant value.
Figure 9. The maximum value of the growth rate ( σ m a x ) on the dispersion relation as a function of the Schmidt number ( S c ) for different Rayleigh numbers ( R a ). In the case of negative Rayleigh numbers, σ m a x continues to increase for large values of S c . While for positive Rayleigh numbers, they approach a constant value.
Symmetry 16 00269 g009
Figure 10. Stream function near a nonaxisymmetric reaction front showing two convective counterrotating rolls. The roll underneath the front is stronger than the one above. The black line represents the position of the front, set at a level where the concentration c is 0.2.
Figure 10. Stream function near a nonaxisymmetric reaction front showing two convective counterrotating rolls. The roll underneath the front is stronger than the one above. The black line represents the position of the front, set at a level where the concentration c is 0.2.
Symmetry 16 00269 g010
Figure 11. Reaction front profiles formed above the critical width. The Rayleigh number is set to R a = 1.2 . The color map represents the dimensionless concentration c. (a) A nonaxisymmetric front develops for a domain width of L x = 30 . (b) An axisymmetric front develops for a domain width of L x = 60 .
Figure 11. Reaction front profiles formed above the critical width. The Rayleigh number is set to R a = 1.2 . The color map represents the dimensionless concentration c. (a) A nonaxisymmetric front develops for a domain width of L x = 30 . (b) An axisymmetric front develops for a domain width of L x = 60 .
Symmetry 16 00269 g011
Figure 12. Front speed for different domain widths ( L x ). There is a region of bistability between the axisymmetric and nonaxisymmetric fronts. The Rayleigh number is set to R a = 1.2 .
Figure 12. Front speed for different domain widths ( L x ). There is a region of bistability between the axisymmetric and nonaxisymmetric fronts. The Rayleigh number is set to R a = 1.2 .
Symmetry 16 00269 g012
Figure 13. Maximum absolute value of the stream function | ψ | m a x for different widths L x of the domain. There is a region of bistability between the axisymmetric and nonaxisymmetric fronts. The Rayleigh number is set to R a = 1.2 .
Figure 13. Maximum absolute value of the stream function | ψ | m a x for different widths L x of the domain. There is a region of bistability between the axisymmetric and nonaxisymmetric fronts. The Rayleigh number is set to R a = 1.2 .
Symmetry 16 00269 g013
Figure 14. Front speed for different domain widths ( L x ). The Rayleigh number is set to R a = 1.2 .
Figure 14. Front speed for different domain widths ( L x ). The Rayleigh number is set to R a = 1.2 .
Symmetry 16 00269 g014
Disclaimer/Publisher’s Note: The statements, opinions and data contained in all publications are solely those of the individual author(s) and contributor(s) and not of MDPI and/or the editor(s). MDPI and/or the editor(s) disclaim responsibility for any injury to people or property resulting from any ideas, methods, instructions or products referred to in the content.

Share and Cite

MDPI and ACS Style

Vilela, P.M.; Guzman, R.; Vasquez, D.A. Thermally Driven Convection Generated by Reaction Fronts in Viscous Fluids. Symmetry 2024, 16, 269. https://doi.org/10.3390/sym16030269

AMA Style

Vilela PM, Guzman R, Vasquez DA. Thermally Driven Convection Generated by Reaction Fronts in Viscous Fluids. Symmetry. 2024; 16(3):269. https://doi.org/10.3390/sym16030269

Chicago/Turabian Style

Vilela, Pablo M., Roberto Guzman, and Desiderio A. Vasquez. 2024. "Thermally Driven Convection Generated by Reaction Fronts in Viscous Fluids" Symmetry 16, no. 3: 269. https://doi.org/10.3390/sym16030269

APA Style

Vilela, P. M., Guzman, R., & Vasquez, D. A. (2024). Thermally Driven Convection Generated by Reaction Fronts in Viscous Fluids. Symmetry, 16(3), 269. https://doi.org/10.3390/sym16030269

Note that from the first issue of 2016, this journal uses article numbers instead of page numbers. See further details here.

Article Metrics

Back to TopTop