Adaptive Finite Element Simulation of Double-Diffusive Convection

Abstract

Double-diffusive convection plays an important role in many physical phenomena of practical importance. However, the numerical simulation of these phenomena is challenging since fine meshes are often required to capture the flow physics. Hence, several different numerical methods have been employed in the past. This work reports the development and application of an adaptive finite element method for the simulation of these phenomena, thereby avoiding the need for the use of very fine meshes over the whole domain. The weak formulation of the conservation equations for mass, momentum, energy and species concentration is used. The Boussinesq approximation relates the density of the fluid to the temperature and/or the species concentration. A second-order backward difference method is used for time discretization and the Galerkin method is employed for spatial discretization. Both adaptive time step and grid refinement techniques are employed, and the code is parallelized using MPI. Three different stabilization methods of the convective-diffusion equations are compared; namely, the streamline upwind Petrov–Galerkin (SUPG) method, and two modified methods aimed at diminishing spurious oscillations that include an artificial diffusion term. This diffusion term may be either isotropic or orthogonal to the streamlines. The addition of artificial isotropic diffusion to the SUPG method provides enhanced stability. The method is applied to double-diffusive finger convection in a sucrose-salt aqueous mixture and a stratified salt solution heated from below. The method accurately reproduces the experimentally observed temporal evolution of the salt fingers in the former case and the location of the interfaces between convective and non-convective zones in the latter.

1. Introduction

Double-diffusive convection was originally described in a paper by Stommel et al. [1], where it was referred to as an oceanographical curiosity. However, only a few years later it was recognized as an important phenomenon with applications in several branches of science. Several early reviews of this topic were published [2,3,4] and books on the subject are presently available [5,6]. Double-diffusive convection is the flow originated by buoyancy in a fluid where two different density gradients are present, with different diffusion rates. These flows occur in a wide variety of fields, such as oceanography, astrophysics, geophysics, solar ponds and crystal growth, among others.

An example of this phenomenon occurs in ocean water with temperature and salinity gradients, e.g., when surface water is warmer and saltier than deeper water. Heat and salt diffuse at different rates—the molecular diffusion coefficient of salt being two orders of magnitude smaller than the thermal diffusion coefficient—yielding double-diffusion convection. Differences in temperature and salinity may cause instabilities, referred to as double-diffusive instabilities, even when the density increases with the depth.

In the case of salty water, in order to achieve stability, density must increase with depth. However, this is not a sufficient condition, since double diffusion phenomena may cause instability, depending on the temperature and salinity gradients. Indeed, four different configurations may be found. When the temperature and the salinity increase with the depth, the temperature gradient has a destabilizing effect, in contrast with the salinity gradient. In this case, a diffusive regime is present, but it may become a convective regime if the potential energy associated with the temperature gradient is high enough. In the opposite case, i.e., when both the temperature and the salinity decrease with the depth, the salinity gradient has a destabilizing effect, in contrast with the temperature gradient, yielding so-called salt fingers. If the temperature increases but the salt concentration decreases with depth, both gradients have a destabilizing effect, leading to a steady convective regime. Finally, if the temperature decreases while the salinity increases with the depth, the system is stable, and only thermal and mass diffusion occurs.

Double-diffusion convection in rectangular enclosures has been investigated by several researchers experimentally, theoretically and computationally. One of the first works on the numerical simulation of double-diffusive convection was reported by Heinrich [7] who applied the finite element method to study a two-dimensional enclosure of salty water. The top and bottom walls were adiabatic, and the vertical walls were maintained at different temperatures. Kazmierczak and Poulikakos numerically investigated transient two-dimensional double diffusion in a horizontal fluid layer with a salt gradient and heated from below [8]. They used a finite difference method and studied the effects of the aspect ratio of the enclosure and the ratio of the thermal to salt Rayleigh numbers. A review of earlier work on the modelling of salt-finger convection was reported by Yoshida and Nagashima [9]. Most of this work was concerned with the finger structures and their temporal variation, and with the variation of the buoyancy–flux ratio as a function of the density–stability ratio. Since then, a significant number of papers have been published on this subject, most of them for salty water with temperature and salinity gradients. Several different numerical methods have been used in the simulations.

Sreenivas et al. [10] used the finite volume method and the SIMPLER algorithm to investigate double-diffusive finger convection in a 2D domain with a step change in salinity and temperature. Simulations were performed for a wide range of buoyancy ratios and thermal Rayleigh numbers, and new insight into the relationships between salt fingers’ width, velocity and fluxes were reported. In a follow-up of this study [11], the morphological evolution of convective structures and salt fingers from low to high thermal Rayleigh numbers was investigated. The authors concluded that the temperature and salinity Rayleigh numbers play an important role in the evolution of the salt fingers. Further research on the influence of the Rayleigh numbers is reported in Ref. [12]. A Hele-Shaw cell [13] was used to experimentally investigate double-diffusive convection. Numerical simulations, based on the same mathematical model and numerical methods used in Refs. [10,11], were also performed. An expression for the growth rate of the salt fingers was derived from the linearized governing equations, which yielded reasonably good predictions.

Gerstenberger et al. [14] used a discontinuous Galerkin method to model viscous fingering instabilities in porous media. The method was applied to the simulation of the dissolution of CO₂ into brine aquifers for sequestration purposes. The same problem was investigated by Shahraeeni et al. [15] who used a combination of a mixed-hybrid finite element method for the fluid velocity and a discontinuous Galerkin method for the species conservation equation. They concluded that 3D fingering can be captured by their method while the finite difference method may be very time-consuming. Norouzi and Shoghi [16] used a spectral method along with the fourth-order Adams–Bashforth technique to solve the equations that describe finger instability in anisotropic porous media.

Seta et al. [17] studied the onset and development of double diffusion instabilities in a ternary mixture. Analytical solutions for pure diffusion and 3D numerical simulations were carried out and validated against experimental data that were obtained using the sliding symmetric tubes technique. The numerical results were obtained using the finite volume method and the OpenFOAM v1812 software. Zhang et al. [18] modelled double diffusion through a fluid-saturated porous interface between a fluid layer overlying a porous layer. They used a single-domain approach, firstly solving the equation for the stream function in the whole domain, before the equations for the vorticity, temperature and salinity were solved in the fluid and porous layers at every time step. Double diffusion instability and salt fingers in porous media were recently investigated by the same authors [19]. They used a novel pore-scale model and solved the governing equations using the lattice Boltzmann method with a multi-relaxation time step. An iterative source-correction immersed boundary method was employed to describe the solid matrix of the porous medium. The role of pores in the evolution of the salt fingers was investigated.

Abdul Hamid and Muggeridge [20] used a second-order finite difference method to investigate the development of fingers over time in unstable, miscible displacement, homogeneous systems. The analysis extended over the lifetime of the fingers, encompassing four regimes, namely their early formation, growth, nonlinear interactions and decay. This work was extended by Kampitsis et al. [21] who developed a high-order conservative dynamic adaptive mesh optimization method to simulate immiscible viscous fingering in porous media. The adaptive mesh was implemented in a double control volume finite element method. They successfully captured the development and growth of the fingers at a much lower computational cost than that required by a method relying on a fixed mesh.

Recently, Ouzani et al. [22] used a finite volume method to numerically investigate double-diffusive salt-fingering convection. A fifth-order scheme was used to discretize the convective terms, a fourth-order centered scheme for the diffusive terms and a third-order TVD Runge–Kutta method for the unsteady term. They showed that the thermal Rayleigh number and the buoyancy ratio have a significant effect on the mixing properties of double-diffusive convection. Liu et al. [23] simulated salt-finger convection in a fluid confined between two parallel and infinitely long plates. They carried out a bifurcation analysis using single-mode equations derived from a truncated Fourier expansion in the horizontal direction. Keable et al. [24] investigated the influence of the viscosity ratio and Peclet number on finger convection in a Hele-Shaw cell. They performed both experiments and computations. A finite difference method was used. They pointed out that a very high mesh resolution is required to simulate the physical diffusion and dispersion and to mitigate numerical diffusion. Ishikawa et al. [25] investigated the transient behavior of double-diffusive convection in two horizontal miscible fluid layers. They used the stream function-vorticity formulation of the governing equations, which were solved using a Fourier spectral method in the horizontal direction and a finite difference method in the vertical direction.

The literature survey shows that several different numerical methods have been employed to model convection in double-diffusive systems and that fine meshes are needed to resolve the small scales. In this paper, an adaptive finite element method is used and fine meshes are restricted to the regions of strong gradients. Hence, the contribution of this work to the literature lies in the development and validation of a finite element method with adaptive time step and mesh refinement, using an open-source finite element library, that allows the simulation of double-diffusive convection, thereby avoiding the need to use very fine meshes extending over the whole computational domain. The performance of two stabilization methods to mitigate undershoots and overshoots of the numerical solution in the presence of strong gradients is compared. Two different problems are addressed: finger convection in a sucrose-salt aqueous mixture (case 1); and double-diffusive convection in a stratified salt solution heated from below (case 2).

2. Mathematical Modelling

2.1. Governing Equations

Double-diffusive convection is governed by the conservation equations of mass, momentum, energy and species mass fraction, which may be written as follows [26]:

$$\frac{\partial \rho }{\partial t}+\nabla \cdot \left(\rho v\right)=0$$

(1)

$$\rho \frac{\partial v}{\partial t}+\rho v\cdot \nabla v=-\nabla p+\nabla \cdot \tau +\rho g$$

(2)

$$\rho \frac{\partial e}{\partial t}+\rho v\cdot \nabla e=-\nabla \cdot q-\nabla \cdot \left(pv\right)-\nabla \cdot \left(\tau \cdot v\right)+\rho g\cdot v$$

(3)

$$\rho \frac{\partial y_i}{\partial t}+\rho v\cdot \nabla y_i=\nabla \cdot \left(\rho D_i\nabla y_i\right)$$

(4)

where the viscous stress tensor, τ, is defined as follows for a Newtonian fluid:

$$\tau =\mu \left(\nabla v+{\left(\nabla v\right)}^T\right)-\frac{2}{3}\mu \left(\nabla \cdot v\right)I$$

(5)

In these equations, v is the velocity vector, e the specific total energy (which is given by the sum of the internal and kinetic energy), y_i the mass fraction of species i, ρ the density, p the pressure, t the time, g the acceleration of gravity vector, q the heat flux vector, I the unit tensor, μ the dynamic viscosity and D_i the diffusion coefficient of species i. Equation (3) is only solved for case 2. Equation (4), which relies on the assumption that the mass diffusion flux is governed by Fick’s law, is solved for the mass fractions of sucrose and salt in case 1, and only for the mass fraction of salt (also referred to as salinity and hereafter denoted by S) in case 2.

The Boussinesq approximation, which assumes that the density is constant except in the buoyancy term of the momentum equations, is used. The density in the buoyancy term is expressed in terms of a Taylor series truncated after the first-order derivative terms, yielding

$$\rho ={\rho }_o+\frac{\partial \rho }{\partial y_1}\left(y_1-y_{1,o}\right)+\frac{\partial \rho }{\partial y_2}\left(y_2-y_{2,o}\right)$$

(6)

for case 1 and

$$\rho ={\rho }_o+{\frac{\partial \rho }{\partial T}|}_p\left(T-T_o\right)+{\frac{\partial \rho }{\partial S}|}_p\left(S-S_o\right)$$

(7)

for case 2. In these equations, ρ_o is the density calculated at reference values of the mass fractions of the species, y_1,o and y_2,o (case 1) and of the temperature, T_o, and salinity, S_o (case 2). The derivatives may be related to the concentration expansion coefficient of the corresponding species, β_Si, and to the thermal expansion coefficient, β_T, as follows:

$${\beta }_{Si}=\frac{1}{\rho }\frac{\partial \rho }{\partial y_i}$$

(8)

$${\beta }_T=-\frac{1}{\rho }{\frac{\partial \rho }{\partial T}|}_p$$

(9)

We further assume that the work due to pressure, viscous and gravitational forces is negligible, that the heat flux is only due to heat conduction and satisfies Fourier’s law, and that the physical properties of the species are constant (apart from the density in the buoyancy term, as stated above). Accordingly, the governing equations may be simplified as follows:

$$\nabla \cdot v=0$$

(10)

$$\frac{\partial v}{\partial t}+v\cdot \nabla v=-\frac{\nabla p}{{\rho }_o}+\nu \Delta v+f$$

(11)

$$\frac{\partial T}{\partial t}+v\cdot \nabla T=\alpha \Delta T$$

(12)

$$\frac{\partial y_i}{\partial t}+v\cdot \nabla y_i=D_i\Delta y_i$$

(13)

where α stands for the thermal diffusivity, D_i for the mass diffusivity of species i and f is defined as

$$f=\left[1+{\beta }_1\left(y_1-y_{1,o}\right)+{\beta }_2\left(y_2-y_{2,o}\right)\right]g$$

(14)

for case 1 and

$$f=\left[1-{\beta }_T\left(T-T_o\right)+{\beta }_S\left(S-S_o\right)\right]g$$

(15)

for case 2.

2.2. Finite Element Discretization

According to the finite element method, both sides of Equations (10)–(13) are multiplied by test functions and all the terms are integrated over the computational domain. Next, the pressure term in the momentum equation and the diffusive terms in the momentum, energy and species equations are integrated by parts. Application of the Gauss divergence theorem yields the weak form of the governing equations. The velocity, pressure, temperature and species mass fractions are approximated by a linear combination of the basis functions and, according to the Galerkin method, the test functions are also expressed as a linear combination of the same basis functions. Bilinear elements were used for the pressure, and biquadratic elements were considered for the velocity, temperature and species mass fractions. The time discretization was carried out using the second-order backward difference method. The integration of the various terms of the governing equations was performed element by element, yielding the local mass, stiffness and load matrices, and then these local matrices were assembled to obtain the global ones. Details may be found in [27].

2.3. Stabilization of Convection-Dominated Equations

The application of the Galerkin finite element method to transport equations yields spurious (unphysical) oscillations whenever the convection is significant compared with the diffusion. In order to overcome this problem, stabilization methods are usually employed. In this work, three different stabilization methods were used; namely, the streamline upwind Petrov–Galerkin (SUPG) method [28], and two modified methods aimed at diminishing spurious oscillations in layers (SOLD) [29], which add an artificial diffusion term to the transport equation. In the SOLD methods, the artificial diffusion term may be either isotropic (SOLD ISO) [30] or orthogonal (SOLD ORT) [31] to the streamlines. The SUPG includes a stabilization parameter that may be calculated using the generalization to multidimensional problems of the parameter originally proposed by Christie et al. [32] for linear elements. Other stabilization parameters have been used in other studies, and that proposed by Tzeduyar and Osawa [33] was employed in the present work. The SOLD methods also include a stabilization parameter. The parameters proposed by Almeida and Silva [34] for the SOLD ISO method and by John and Knobloch [29] for the SOLD ORT method were used here. The performance of these three stabilization methods (SUPG, SOLD ISO and SOLD ORT) was investigated for the classical problem of pure convection of a scalar discontinuity by a uniform velocity field by Milhazes [35].

2.4. Adaptive Mesh Refinement and Time Discretization

The adaptive mesh refinement was implemented using the Kelly error estimator (KEE) proposed by Kelly et al. [36], which is based on gradient recovery. The KEE estimates the error in each element through the jump in the gradient of a scalar on the interelement boundaries. The adaptive refinement is controlled by three parameters. The first parameter is the percentage of elements to be refined, coarsened and maintained unchanged. The second parameter is the number of time steps between successive adaptive passes. The last parameter is the scalar used to estimate the error. In the present work, the selected scalar is the density, the gradient of which is determined from Equation (6) or (7) according to the Boussinesq approximation:

$$\nabla \rho ={\rho }_o\left({\beta }_1\nabla y_1+{\beta }_2\nabla y_2\right)$$

(16)

$$\nabla \rho ={\rho }_o\left(-{\beta }_T\nabla T+{\beta }_S\nabla S\right)$$

(17)

The time step was adaptively controlled based on the method proposed by Jonhson et al. [37]. At mth time instant, and for a variable ϕ, the time step, Δt^m, should satisfy the following condition:

$$C\mathsf{\Delta }{t}^m{\Vert {\varphi }^m-{\varphi }^{m-1}\Vert }_{\infty }\le \epsilon$$

(18)

where ε is the tolerance, taken as 10⁻⁵, and C is the Courant number defined as

$$C={\Vert u\Vert }_2\frac{\mathsf{\Delta }{t}^m}{h}$$

(19)

Here, ${\Vert u\Vert }_2$ is the maximum velocity in the Euclidian norm and h is a characteristic element size, taken as the diagonal of the smallest element. If the condition given by Equation (18) is satisfied at the mth time instant for all equations, then the time step is increased by a small factor, taken as 1.005, and the algorithm advances to the next time instant. Otherwise, if the condition described in Equation (18) is not fulfilled, the time step is reduced by 5% and the solution is recalculated.

2.5. Solution Algorithm

The pressure and the velocity are coupled in the Navier–Stokes equations. In this work, a projection method was used to solve these equations. These methods handle the Navier–Stokes equations in such a way that a vectorial advection-diffusion and a Poisson equation are solved.

Let us assume that the initial velocity and pressure fields are equal to zero, as well as an auxiliary scalar φ. The second-order backward difference method is considered for time discretization. The projection method [38] involves four steps. In the first step, the velocity and the pressure are extrapolated as follows:

$$v^{*}=2v^m-v^{m-1},p^{*}=p^m+\frac{4}{3}{\phi }^m-\frac{1}{3}{\phi }^{m-1}$$

(20)

In the second step, the velocity field v^m+1 is found by solving the linear equation derived from the skew-symmetric form of the momentum equation [39]:

$$\frac{1}{2\mathsf{\Delta }t}\left(3v^{m+1}-4v^m+v^{m-1}\right)+v^{*}\cdot \nabla v^{m+1}+\frac{1}{2}\left(\nabla \cdot v^{*}\right)=-\frac{\nabla p^{*}}{{\rho }_o}+\nu \mathsf{\Delta }{v}^{m+1}+f^{m+1}$$

(21)

In the third step, the following differential equation is solved:

$$\mathsf{\Delta }{\phi }^{m+1}=\frac{3}{2\mathsf{\Delta }t}\nabla \cdot v^{m+1}$$

(22)

In the last step, the pressure is corrected as follows:

$$p^{m+1}=p^m+{\phi }^{m+1}-\mu \nabla \cdot v^{m+1}$$

(23)

The calculation of the scalar fields (species mass fractions and/or temperature) at time t^m+1 requires knowledge of the velocity at that time. Similarly, the buoyancy term f in the Navier–Stokes equations requires knowledge of the species mass fractions and/or temperature at the same time. The coupling of the velocity, species mass fractions and/or temperature is achieved using a linear extrapolation of these variables at time t^m to estimate the unknown values at time t^m+1, according to the following algorithm (where T*, S* and v* denote extrapolated values):

(i).

Calculate y₁^m+1 and y₂^m+1 (or T^m+1 and S^m+1) using v*

(ii).

Calculate the buoyancy term, f^m+1, using y₁* and y₂* (or T* and S*)

(iii).

Calculate the solution of the Navier–Stokes equations using the projection method and f^m+1

(iv).

t^m+1 = t^m + Δt^m

(v).

Return to step (i)

2.6. Parallelization of the Code

The deal.II open-source finite element library [40] allows the parallelization of the code using threading building blocks (TBB). TBB is a software template library for task-based shared memory parallel programming on multicore processors, which was developed by Intel. The mesh is partitioned into a number of blocks equal to the number of available processors, and every element is assigned to a block. The basic linear algebra operations are parallelized using the available libraries, but the user needs to define other types of parallelization, such as the assembling of the global matrices and the solution of the system of linear algebraic equations. The calculation of the local element matrices and the assembling of these matrices are parallelized.

3. Results and Discussion

The code used in the present simulations was formerly applied to the simulation of a small-scale salt gradient solar pond and the predictions were compared with both experimental data and numerical results obtained using the commercial software Ansys-Fluent, as reported by Giestas et al. [41]. Two different double-diffusive convection problems are addressed in the present work. In the first one, a sucrose-salt aqueous mixture is considered, while in the second one, a salt water solution with a temperature gradient is simulated.

3.1. Case 1: Double-Diffusion Fingers Convection in a Sucrose-Salt Aqueous Mixture

A hydrostatic system may be formed when a low-density fluid layer rests on top of a high-density fluid layer. However, if the fluid has a component whose concentration is greater on the top than on the bottom layer, the different rates of mass diffusion lead to the formation of salt fingers and promote hydrodynamic instability. The same phenomenon is observed in a system of two aqueous solution layers with different molecular diffusion. Several authors studied experimentally the formation of salt fingers in a two-dimensional system by placing the aqueous solutions between two acrylic parallel plates very close to each other and using a dye for visualization purposes. This system is referred to as a Hele-Shaw cell [13]. The flow is identical to a two-dimensional fluid flow through a highly porous medium.

The experiments conducted by Cooper et al. [42] are numerically simulated in this section. The Hele-Shaw cell is 0.17 m long and 0.31 m high. The bottom half contains a salt solution, and the top half contains a sucrose solution. The two aqueous solutions are separated by a thin plastic slider, which is removed at the beginning of the each experiment.

The calculations were carried out by adding a term that simulates the porosity to momentum equation, which is now written as follows:

$$\frac{\partial v}{\partial t}+v\cdot \nabla v=-\frac{\nabla p}{{\rho }_o}-\nu \frac{v}{K}+\nu \mathsf{\Delta }v+f$$

(24)

where K is the intrinsic permeability of the medium, which is taken as b²/12 [42]. Here, b denotes the distance between the two plates of the Hele-Shaw cell. The following data were used in the calculations: ρ_o = 1.05 × 10³ kg/m³, ν = 1.025 × 10⁻⁶ m²/s, D_S = 1.57 × 10⁻⁹ m²/s for salt; D_S = 5.53 × 10⁻¹⁰ m²/s for sucrose; β_S = 0.755 for salt; β_S = 0.415 for sucrose.

Dirichlet homogeneous boundary conditions for the velocity and homogenous Neumann boundary conditions for the mass fraction of the species are prescribed along the boundary. The initial velocity is equal to zero. The initial salinity is equal to 10.88 g/kg in the lower half of the domain and zero elsewhere, while the initial concentration of sucrose is constant in the top half, its value being dependent on the prescribed density stability ratio, and zero elsewhere. The initial mesh has 80 × 144 elements.

The initial salt concentration is 13.81 g/kg, which corresponds to a density stability ratio Rρ = 1.4. This parameter is defined as Rρ = (β_S1 Δc_o,1)/(β_S2 Δc_o,2), where Δc_o denotes the difference in initial concentration between the two layers, and subscripts 1 and 2 denote the faster and slower diffusing species, respectively. The instabilities that originate salt fingers begin at a time that depends on Rρ. In this case, the first instabilities are observed at 460 s. As they grow, a structure of well-spaced and organized salt fingers develops from the interface between the two aqueous solutions, and the fingers become progressively longer and more well-defined. When a tubular structure appears, the fingers interact among them, causing a reorganization of the initial uniform structure and evolving to a less uniform one, that exhibits differences in the fingers’ length.

The experimental and computational results for sucrose concentration, normalized by the maximum sucrose concentration, at t = 52 min, are displayed in Figure 1 and reveal a good qualitative agreement. In both cases, the tubular structures are visible, as is their interaction. A few fingers are longer than others. During the transient evolution of the double-diffusive convection, new small fingers are formed at the central region, referred to as generation zone [42]. The new fingers may grow independently, merge with adjacent fingers or entrain previously formed fingers. At t = 166 min, the length of the fingers is greater than before. A few branches in the structures appear; the interaction between neighbouring structures is more visible and new fingers merge with previously existing ones. At t = 52 min, the average size of the fingers is similar in the simulations and experiments, while at t = 166 min the average length is slightly underestimated below the generation zone. A zoom of the generation zone is shown In Figure 2, which also shows the adaptive mesh. The mesh is more refined in the fingers’ region, as expected since the gradient of the sucrose concentration is greater in that region and coarser elsewhere due to the lower concentration gradient.

The Increase of Rρ causes a delay in the evolution of the structures observed above and these structures are more diffuse. The fingers grow slower, in a direction that is closer to the vertical direction and exhibit diffuse regions between them. The branching formerly observed is more sporadic. In the case of Rρ = 2.8, which corresponds to an initial sucrose concentration of 7.1 g/kg, the first instabilities are observed at about t = 2250 s. The tubular structures grow slowly, and diffusion prevails. Figure 3 shows the predicted and experimentally visualized normalized sucrose concentration at two different times. At t = 522 min, the growth of the structures is more uniform than that observed for the lower Rρ and there is more diffusion between the small fingers, which do not exhibit any branches. The thickness of the generation region is greater, and the length of the fingers is smaller in comparison with that observed in the case of Rρ = 1.4. At t = 2993 min, the growth of the tubular structures reaches a maximum, and the fingers achieve a maximum length of 8 cm in the experiments, and a little less in the computations.

The performance of the three different stabilization methods employed in this study was compared for the case Rρ = 1.4. In addition to the SOLD ISO method used in the calculations shown in Figure 1, Figure 2 and Figure 3, the SUPG and SOLD ORT methods were also used, along with the stabilization parameter proposed by Tezduyar and Osawa [33]. In order to investigate the occurrence of overshoots and undershoots, the absolute value of the difference between the maximum and the minimum sucrose concentration normalized by the corresponding initial difference was evaluated. Values of 4.7%, 0.4% and 0.3% were obtained for SUPG, SOLD ISO and SOL ORT, respectively. In the former case, the fingers are well-defined, but oscillations in the maximum and minimum values of the sucrose concentration were found at the beginning of the instabilities. The SOLD ISO method also yields solutions with well-defined fingers, as shown in Figure 1 and Figure 2, and the overshoots/undershoots are largely suppressed. The numerical solution is physically realistic, but the computational cost is greater than that of the SUPG method. When the SOLD ORT method is employed, the fingers are a little more diffuse, but the solution is qualitatively similar to that of the SOLD ISO method. The overshoots/undershoots are marginal.

Additional calculations were performed using the stabilization parameter of Christie et al. [32] instead of that of Tezduyar and Osawa [33]. In this case, the solutions obtained using the SOLD ISO and SOLD ORT methods do not exhibit fingers, i.e., they are too diffusive, which means that the stabilization parameter from Ref. [32] does not perform well. The SUPG method is less diffusive and allows the formation of fingers, but these are incipient and rather diffuse, causing the end of convection after the beginning of the instabilities. Accordingly, the numerical solution is very poor, as revealed by the results obtained at t = 52 min and t = 166 min (see Figure 4), which are very different from the experimental data and the numerical solution shown in Figure 1. These results are somewhat similar to those of Rρ close to the critical value beyond which only diffusion occurs, which was found to be equal to 3.0 in the experiments of Cooper et al. [42]. Therefore, the stabilization parameter proposed in Ref. [33] should be used instead of that proposed in Ref. [32].

3.2. Case 2: Double-Diffusion Fingers Convection in a Two-Layer Heat-Salt System

The second test case consists of double diffusion in salt-stratified water in a rectangular enclosure whose lower surface is heated. Experiments were formerly performed in a laboratory plexiglass box with 5 cm width and 10 cm height [43]. The initial velocity is zero and the initial temperature is equal to the ambient temperature (22.95 °C). The initial salinity varies linearly from 51.2 g/kg at the bottom surface to 12.3 g/kg at a height of 8.5 cm and is constant from there to the top surface. In this problem, the temperature gradient has a destabilizing effect, but the salinity gradient does not. Therefore, either a diffusive or convective regime may occur, but salt fingers are not formed.

Homogeneous Dirichlet boundary conditions are used for the velocity. The temperature at the bottom boundary is equal to 50 °C, but the application of a Dirichlet boundary condition for the temperature at this boundary leads to an unstable numerical solution due to the discontinuity between the temperature of the boundary and the initial temperature of the fluid. In order to overcome this problem, a Neumann boundary condition is used. In the initial transient period, only conduction takes place, and the heat flux is determined from the classical equation for a semi-infinite medium subjected to a sudden change of temperature at the boundary [44]. After convection begins, the heat flux is computed from Newton’s law of cooling with a convective heat transfer coefficient estimated from an empirical correlation [45]. A Robin boundary condition for the temperature is used at the side and top boundaries, with the heat transfer coefficients estimated again from empirical correlations. The thermal resistance of the wall is neglected. A homogeneous Neumann boundary condition is used for the salinity at the side and top walls. The initial mesh has 70 × 140 elements. The time step size is limited to a maximum of 0.06 s.

Figure 5, Figure 6, Figure 7 and Figure 8 show the results obtained at 580, 1070, 2470 and 4680 s. In these figures, LCZ and NCZ denote the lower convective zone and the non-convective zone, respectively, while CZ1 and CZ2 denote two smaller and temporary convective zones. Only the lower part of the enclosure is shown for clarity purposes, so that the UCZ (upper convective zone) is not visible. The predicted velocity field and the contours of temperature and salinity are shown over Schlieren images. These images, which were experimentally obtained [43], provide information on the variation of the density gradient through the change in the light intensity (greyscale). The darker regions in the Schlieren images correspond to greater changes in density, associated with steep gradients of temperature and/or salinity, and are typically associated with the transition between different flow zones.

At t = 580 s, a disturbance can be observed in Figure 5, which is due to a plume of rising hot fluid that has a high enough momentum to penetrate through the interface between the LCZ, in the lower region of the domain, and the NCZ, in the upper region, and extend about 0.5 cm above it. In this process, the plume transports heat and mass from the LCZ to the NCZ. This phenomenon is accompanied by a variation in density in the immediate surroundings of the disturbance, as can be seen from the Schlieren image. The temperature and the density are approximately uniform in the LCZ and decrease sharply across the interface. The instability mechanism occurs when the LCZ is still relatively small, being unable to absorb the kinetic energy of the hot plume and prevent the breakdown of the LCZ–NCZ interface. Figure 5 also shows the fluid moving down to the bottom of the enclosure in the vicinity of the walls, due to lateral cooling. This mechanism destabilizes the interface throughout its extension, while the salt convection contributes to decreasing the salinity in the LCZ.

The disturbance caused by the rising plume shown In Figure 6 disappears over time as the heat diffuses in the NCZ. Another interface destabilization mechanism is illustrated in Figure 6, at t = 1070 s. In this case, the kinetic energy associated with the salinity transport in the vicinity of the interface is sufficiently high to break it down, allowing the entrainment of fluid immediately above the NCZ into the LCZ. This mechanism locally destabilizes the interface and contributes to decreasing the salinity in the LCZ. The transport of heat and mass along the walls towards the bottom is also visible here. The height of the LCZ is a little greater than that observed at t = 580 s.

The system continues to evolve with the sporadic appearance of disturbances, such as that shown in Figure 5. Figure 7 illustrates the stratification mechanism due to double diffusion with side cooling at t = 2470 s. It leads to two new thin but well-defined convective zones, CZ1 and CZ2, which appear between the LCZ and the NCZ, and whose height is similar, while the height of the LCZ has increased, as expected. Vortices close to the left wall appear in the two new convective zones, as well as vortices next to the right wall. Similar vortices have been observed in experiments of double-diffusive convection with side cooling [46,47]. The temperature decreases in the vicinity of the walls due to heat losses. However, both the temperature and the salinity remain approximately constant in the LCZ, as well as in the CZ1 and CZ2, and decrease sharply at the interfaces between the different zones (LCZ and CZ1; CZ1 and CZ2; and CZ2 and NCZ).

At t = 4680 s, the CZ2 previously observed merged with the NCZ and only the CZ1 with about 1 cm height is still present between the LCZ and the NCZ (see Figure 8). The vortices in that zone move in the clockwise direction close to the right wall and in the anticlockwise direction next to the left wall, as expected in the case of lateral cooling. The salinity is approximately homogeneous in the CZ1, apart from small variations along the vortices. The heat and salinity transport mechanism along the vertical walls towards the bottom is still present.

Figure 9 shows the time evolution of the location of the interfaces between different zones, which are defined by steep gradients in the density in the vertical direction. In this figure, LCZ, CZ1 and CZ2 denote the location of the upper interfaces between these zones and the neighbouring ones, while UCZ denotes the location of the interface between the NCZ and the UCZ. The height of the LCZ grows significantly in the first few minutes and then continues to increase more slowly up to the end of the simulation. The theory shows that in the case of insulated side walls, this height is proportional to the square root of the time [2]. However, in the present case there is heat loss through the side boundaries and therefore a lower growth rate is expected. The predicted height of the LCZ was fitted by a power function of the time, as shown in Figure 9, and it was found that the exponent of the fitting function is 0.42, which is consistent with the expected magnitude of the temporal evolution of that height. Shortly after 40 min, the first instabilities due to side cooling appear, which lead to the formation of a new convective zone, CZ1. At about 52 s, another convective zone, CZ2, is formed. However, in the period ranging from 70 to 92 min, these two convective zones merge into a single one. Then, that convective zone breaks down, leading again to two convective zones that prevail until about 104 min. Thereafter, CZ2 disappears and only CZ1 persists until it also disappears at about 150 min. The location of the interface of the UCZ remains approximately constant for about 80 min and then slowly moves down.

The temporal evolution of the temperature and salinity at different depths is plotted in Figure 10 where z stands for the depth, i.e., the distance from the top surface to the horizontal plane under consideration. The temperature and salinity at a given z represent the average temperature and salinity, respectively, at that z location and at the time under consideration. The temperature increases at depths greater than or equal to 4 cm until a peak is attained. At a later time, the temperature decreases, and its evolution is identical to that observed at greater depths, i.e., the temperature is uniform at depths greater than that under consideration, which corresponds to the LCZ. In fact, the temperature in the LCZ is approximately uniform, except in the close vicinity of the walls, as illustrated in Figure 5, Figure 6, Figure 7 and Figure 8. At a depth of 3 cm, the evolution is similar to that observed at higher depths, except that there is a decrease between 4.5 and 5 h. At lower depths, the temperatures tend towards that observed at the top surface because the temperature is approximately uniform in the UCZ, except in the region close to the walls.

The evolution of the salinity at different depths, also shown in Figure 10, exhibits features similar to those observed for the temperature. The main difference is that the salinity at a given depth initially remains approximately constant and after some time increases or decreases suddenly, depending on the depth, and then follows an evolution similar to that at the greatest or lowest depths, corresponding to the LCZ and UCZ, respectively. All the lines converge to the same salinity of 2.85% at about 9.5 h, which is in agreement with the value reported in Tavares [43] for the homogenization of the aqueous solution.

4. Conclusions

Double-diffusive convection was numerically simulated in a two-dimensional rectangular domain using an adaptive finite element method to solve the governing equations for mass, momentum, energy and species conservation. Several stabilization methods were investigated. The SUPG method is prone to unphysical overshoots and undershoots, while the SOLD ISO and SOLD ORT methods mitigate these effects, but at the expense of the introduction of artificial diffusion. The latter methods are preferable in modelling transport phenomena, particularly double-diffusive convection, since the buoyancy term in the momentum equations has an opposite influence on the temperature and salinity. The stabilization parameter also plays a key role. Two different configurations were considered: fingers convection in a sucrose-salt aqueous solution, and double diffusion in a heated salt solution. In the first case, the formation and evolution of salt fingers in a Hele-Shaw cell were satisfactorily predicted in comparison with available experimental data. In the second case, the locations of the interfaces between different convective and diffusive zones were compared with Schlieren images and good qualitative agreement was observed. Overall, the unsteady heat and mass transfer phenomena under investigation were successfully simulated using a finite element method with adaptive mesh and time step techniques.