Soumaya Ouadhani

Ali Abdennadher

Abdelkader Mojtabi

Analytical and numerical stability analysis of Soret-driven convection in a horizontal porous layer: The effect of vertical vibrations

The authors studied the effect of vertical high-frequency and small-amplitude vibrations on the separation of a binary mixture saturating a porous cavity. The horizontal bottom plate was submitted to constant uniform heat flux and the top one was maintained at constant temperature while no mass flux was imposed. The influence of vertical vibrations on the onset of convection and on the stability of the unicellular flow was investigated for positive separation ratio ψ. The case of high-frequency and small-amplitude vibrations was considered so that a formulation using time averaged equations could be used. For an infinite horizontal porous layer, the equilibrium solution was found to lose its stability via a stationary bifurcation leading to unicellular flow or multicellular one depending on the value of ψ. The analytical solution of the unicellular flow was obtained and separation was expressed in terms of Lewis number, separation ratio and thermal Rayleigh number. The direct numerical simulations using the averaged governing equations and analytical stability analysis showed that the unicellular flow loses its stability via oscillatory bifurcation. The vibrations decrease the value of ψ uni, which allows species separation for a wide variety of binary mixtures. The vibrations can be used to maintain the unicellular flow and allow species separation over a wider range of Rayleigh numbers. The results of direct numerical simulations and analytical model are in good agreement.

Introduction

Double-diffusive convection caused by temperature and concentration gradients in a porous medium is a classical example of the problems that reveal the interaction of different instability mechanisms. Reviews of recent developments and publications in this field can be found in [START_REF] Nield | Convection in Porous Media[END_REF]. A great number of works in this field of research are devoted to investigation into high-frequency and small amplitude vibration. Khallouf et al. [START_REF] Khallouf | [END_REF] considered a square differentially heated cavity filled with a porous medium saturated by a pure fluid and subjected to linear harmonic oscillations in the vertical direction. In their study, the authors used a Darcy-Boussinesq model and a direct formulation. Gershuni et al. [3] have studied the stability of a substantially linear mechanical equilibrium solution for an horizontal layer of a binary mixture subjected to ver-Contribution to the Topical Issue "Non-isothermal transport in complex fluids", edited by Rafael Delgado-Buscalioni, Mohamed Khayet, José María Ortiz de Zárate and Fabrizio Croccolo.

a e-mail: mojtabi@imft.fr tical thermal gradient under the effect of high-frequency vibrations. The thermovibrational convection in a fluid received more attention than the thermovibrational convection in a porous medium. A good summary of this work was achieved by Gershuni et al. [START_REF] Gershuni | Thermal Vibrational Convection[END_REF][START_REF] Gershuni | [END_REF].

In the case of a horizontal porous layer saturated by a pure fluid, heated from below or from above Zenkovskaya et al. [6] and Bardan et al. [7] used the Darcy model including the non-stationary term and adopted the timeaveraged equations formulation to study the influence of high-frequency and small-amplitude vibrations on the onset of convection. They found that vertical vibrations stabilize the rest solution. Charrier-Mojtabi et al. [8] investigated the influence of vibrations on Soret-driven convection in a horizontal porous cell heated from below or from above. They showed that the vertical vibrations had a stabilizing effect while the horizontal vibrations had a destabilizing effect. Elhajjar et al. [9] have studied the influence of the vertical vibration of high frequency and small amplitude on separation for a cell with horizontal walls maintained at constant but different temperatures. The linear stability of the unicellular flow obtained was analyzed analytically and numerically. They concluded that the unicellular flow associated with a field of laminate concentration leads to the horizontal separation of binary mixture components. They found that the vibrations have a stabilizing effect, leading to an increase of the critical value of the Rayleigh number corresponding to the transition between the unicellular and multicellular flow. Consequently vertical vibration allows species separation for a wide variety of binary mixtures. The validity of the average approach for the description of vibrational convection in porous medium differs from the case with homogeneous fluid and is not widely known. This point is developed in more details in the paper by Lyubimov et al. [10]. The authors investigated the onset and nonlinear regimes of convection in a two-layer system composed of a horizontal pure fluid layer and a fluid-saturated porous layer subjected to the gravity field and high-frequency vibrations.

In the present paper, we use the same formulation as the one used by Charrier-Mojtabi et al. [8] for a shallow porous cavity saturated by a binary mixture. The cavity bounded by horizontal infinite or finite boundaries. The bottom horizontal plate is submitted to constant uniform heat flux and the top one is maintained at constant temperature, while no mass flux is imposed. We verify that it is possible to carry out the species separation of a binary mixture in this geometrical configuration. We consider the case of high-frequency and small-amplitude vibrations, so that a formulation using time averaged equations can be employed. The results of the linear stability analysis of the mechanical equilibrium and the unicellular flow in an infinite porous layer is studied analytically and using spectral method.

Mathematical formulation

We consider a rectangular cavity with aspect ratio A = L/H, where H is the height of the cavity along the vertical axis and L is the width along the horizontal axis (fig. 1). The aspect ratio is assumed infinite in the stability analysis. The cavity is filled with a porous medium saturated by a binary fluid for which the Soret effect is taken into account. Dirichlet and Neumann boundary conditions for temperature are applied to the horizontal walls (z = 0, z = H). The vertical walls (x = 0, x = L) are impermeable and adiabatic. All the boundaries are assumed rigid. The cavity is subjected to linear harmonic oscillations in the vertical direction (amplitude b and dimensional frequency ). For the governing equations, we adopt the Boussinesq approximation and Darcy equation for which the non-stationary term is taken into account.

We set all the properties of the binary fluid constant except the density ρ in the buoyancy term, which depends linearly on the local temperature and mass fraction:

ρ = ρ r [1 -β T (T -T r ) -β C (C -C r )], (1)
where ρ r is the fluid mixture density at temperature T r and mass fraction C r . β T and β C are the thermal and concentration expansion coefficients respectively. When we consider the referential related to the oscillating system, the gravitational field g is replaced by g + b 2 sin( t ) e z where e z is the unit vector along the vertical axis (vibration axis) and t the dimensional time.

The reference scales are H for the length, σH 2 /a for the time, with a = λ * /(ρc) f and σ = (ρc) * /(ρc) f (where λ * and (ρc) * are, respectively, the effective thermal conductivity and volumetric heat capacity of the porous medium), a/H for the velocity (a is the effective thermal diffusivity), ΔT = q H/λ * for the temperature and

ΔC = -ΔT C i (1 -C i ) D *
T /D * for the mass fraction, where C i , D * T , D * are the initial mass fraction, the thermodiffusion and the mass diffusion coefficients of the denser component, respectively. The dimensionless temperature and mass fraction are, respectively, defined by

T = (T -T 2 )/ΔT , C = (C -C i )/ΔC.
Thus, the dimensionless governing conservation equations for mass, momentum, energy and chemical species, where the Soret effect is taken into account, are

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∇ • V = 0, B∂ t V + V = -∇P + Ra(T + ψC)(1 -R sin(ωt)) e z , ∂ t T + V • ∇T = ∇ 2 T, ∂ t C + V • ∇C = 1 Le (∇ 2 C -∇ 2 T ), (2)
where B = Da(ρc) f /(ρc) * Pr is the inverse Vadasz number (B = 1/V a), Da = K/H 2 is the Darcy number and K the permeability of the porous medium.

The problem under consideration depends on eight non-dimensional parameters: the thermal Rayleigh number,

Ra = Kgβ T HΔT /aν, R = bω 2 /g, the separation ra- tio ψ = -(β c /β T ) (D T /D) C i (1 -C i ), the Lewis number Le = a * /D * , the normalized porosity = * (ρc) f /(ρc) *
(where * is the porosity), the dimensionless frequency ω = ω 2 (ρc) * H 2 /a * (ρc) f , the aspect ratio A and the factor B. The dimensionless boundary conditions are

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ V • n = 0, ∀M ∈ ∂Ω, ∂C ∂z = ∂T ∂z = -1, z = 0 ∀x ∈ [0, A], ∂C ∂z = ∂T ∂z , z= 1 ∀x ∈ [0, A], T = 0, z= 1 ∀x ∈ [0, A], ∂C ∂x = ∂T ∂x = 0, x= 0, A ∀z ∈ [0, 1]. (3)
In the momentum equation the term B∂V /∂t is usually neglected since B is of order 10 -6 . However, in our problem, high-frequency vibrations cause very large accelerations, making it necessary to consider this non-stationary term [6].

The averaged equations

In the limiting case of high-frequency and small-amplitude vibrations, the averaging method can be applied to study thermal vibrational convection. According to this method, each field (V , P , T , C) is subdivided into two parts: the first part varies slowly with time (i.e., the characteristic time is large with respect to the period of the vibrations) and the second one varies quickly with time (i.e., the characteristic time is of the order of magnitude of the vibrational period):

(V , P, T, C) = (V * , P * , T * , C * )(t) + (u , p , θ , c )(ω, t).
(4) Here, V * , P * , T * , C * are the averaged fields (i.e., the mean value of the field calculated over the period τ = 2π/ω) of the velocity, pressure, temperature and mass fraction. The decoupling between the pulsational parts of the velocity and the pressure is obtained by using a Helmholtz decomposition:

(T * + ψC * )e z = W + ∇ξ, ( 5
)
where W is the solenoidal part of (T * + ψC * )e z satisfying ∇ • W = 0. Thus, the averaged flow equations are

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ∇ • V * = 0, B∂ t V * + V * = -∇P * + (Ra(T * + ΨC * ) +Rv ∇T * + Ψ * ∇C * • W )e z , ∂ t T * + (V * • ∇)T * = ∇ 2 T * , ∂ t C * + (V * • ∇)C * = 1 Le (∇ 2 C * -∇ 2 T * ), (T * + ψC * )e z = W + ∇ξ, ∇ • W = 0. (6)
In addition to the boundary conditions ( 6)-( 10) applied to the mean fields, we assume that

W • n = 0 on ∂Ω. ( 7
)
The modified vibrational Rayleigh number Rv =

Ra 2 R 2 B 2(B 2 ω 2 +1
) characterizes the intensity of the vibrations.

Linear stability of the equilibrium solution in an infinite horizontal porous layer

The stability of the equilibrium solution was studied by Charrier-Mojtabi et al. [8]. They restricted their study to the case Le = 2 for which the fluid considered is in the gaseous state, and so the Dufour effect should be taken into account. We extended this study to the case of a high Lewis number and we focused on the transition from the equilibrium solution to the unicellular flow obtained for binary mixtures. This problem admits a mechanical equilibrium solution characterized by

V * 0 = 0, T * 0 = -z + cst, C * 0 = -z + 1 2 , W 0 = 0. (8)
In order to analyze the stability of this conductive solution, we introduce the velocity stream function ψ , temperature θ, the transformation η = cθ, and finally, the solenoidal stream function ϕ , characterized by v x = ∂ψ ∂z , v z = -∂ψ ∂x and w x = ∂ϕ ∂z , w z = -∂ϕ ∂x . We assumed that the perturbations ( ψ (z), θ(z), c(z), ϕ (z)) are small. The perturbation quantities are chosen as follows:

(ψ , θ, c, ϕ ) = ψ (z), θ(z), c(z), ϕ (z) e Ikxx+σt . ( 9
)
We introduce a new function η = cθ. The system of equations can be written as

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (Bσ + 1)(D 2 -k 2 ) ψ + Raik((1 + ψ) θ + ψ η) -Rv ϕ k 2 1 + ψ = 0, (D 2 -k 2 ) θ -σ θ -ik ψ = 0, (D 2 -k 2 ) η/Le -σ(η + θ) -ik ψ = 0, (D 2 -k 2 ) ϕ + ik((1 + ψ) θ + ψ η) = 0, (10)
where D = ∂ ∂z , k is the wave number in the horizontal (ox) direction, I 2 = -1, and σ is the temporal exponential growth rate of perturbation. The corresponding boundary conditions are

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ z = 0: ψ = ϕ = ∂ θ ∂z = ∂ η ∂z = 0; z = 1: ψ = ϕ = θ = ∂ η ∂z = 0. ( 11
)

Limiting case of the long-wave mode instability

We study the special case of the long-wave mode theoretically. In some related studies in fluid media, Razi et al. [11] showed that the asymptotic analysis results in a closed form relation for the stability threshold. In order to obtain such a relation, a regular perturbation method with the wave number as the small parameter is performed (for simplifying the procedure), we drop the tilde symbol and we expand (ψ , θ, η, ϕ , σ) as

(ψ , θ, η, ϕ , σ) = ∞ n=0 k n (ψ n , θ n , η n , ϕ n , σ n ), (12)
with boundaries conditions:

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ z = 0: ψ n = ϕ n = ∂θ n ∂z = ∂η n ∂z = 0; z = 1: ψ n = ϕ n = θ n = ∂η n ∂z = 0. (13)
By substituting expressions in the amplitude equations resulting from the linear stability analysis and factoring the same order of k, we find a sequential system of equations: For the zeroth order (k 0 ):

ψ 0 = 0, ϕ 0 = 0, θ 0 = 0, η 0 = cst, σ 0 = 0. ( 14
)
At order 1, we get:

θ 1 = 0, ψ 1 = -IRaψη 0 (z 2 -z)/2, η 1 = cst, ϕ 1 = -Iψη 0 (z 2 -z)/2, σ 1 = 0 (15)
and at order 2:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ θ 2 = Raψη 0 (z 4 -2z 3 + 1)/24, ψ 2 = -IRaψη 1 (z 2 -z)/2, η 2 = η 0 z 2 2 (Le σ 2 + 1) + RaLeψη 0 24 z 3 (z -2), ϕ 2 = -Iψη 1 (z 2 -z)/2. ( 16
)
After involving the solvability condition, we find:

σ 2 = 1 Le -ψ Ra 12 . ( 17
)
We note that σ 2 is real indicating that the conductive solution loses stability through a stationary bifurcation.

For the marginal stability σ 2 is set equal to zero and we obtain the Rayleigh number:

Ra = 12 ψLe . ( 18
)

Stability analysis results for arbitrary values of wave number

The disturbances are developed in terms of polynomial functions verifying all the boundary conditions except those along the inner plates: The Tau spectral method used is similar to the Galerkin method. The test functions used for ψ , θ, ϕ and η verify all the boundary conditions.

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ψ = N n=1 a n (1 -z)z n , θ = b 1 z 2 - 2 3 z 3 + b 2 z 3 -1 + N n=1 b n+2 (1-z) 2 z n+1 , η = d 1 + d 2 z 2 - 2 3 z 3 + N n=1 d n+2 (1 -z) 2 z n+1 , ϕ = N n=1 e n (1 -z)z n . ( 19
)
For ψ > 0 the first bifurcation is stationary. The results are presented in the stability diagram Ra c = f (ψ) and kc = f (ψ). Figures 2 and3 illustrate the effect of vibration on the onset of convection. For the case Le = 100, = 0.5, B = 10 -6 and for Rv = 0, 10, 50, it can be noted that Ra c increases with Rv whereas k c decreases with Rv. Figure 4 shows the influence of vibrations on ψ uni , the separation ration beyond which the critical wave number vanishes (k c = 0), decreases. We note that ψ uni decreases and becomes close to zero when Rv increases. So by adding vibrations, we can separate most binary mixtures. In the case Rv = 0, we obtain ψ uni =

Analytical solution of the unicellular flow

In the case of a shallow cavity (A 1) in order to solve the problem analytically, the parallel flow approximation is considered [9,12]. The streamlines are assumed to be parallel to the horizontal walls except for the vicinity of the vertical walls. In this case, the vertical component of velocity can be neglected. The temperature and mass fraction profiles are written as the sum of two terms: the first one defining the linear longitudinal variation and the second one giving the transverse distribution: The basic flow, denoted with a subscript 0, is then given as follows:

V 0 = U (z)e x , T 0 = bx + h(z), C 0 = mx + g(z), W 0 = W 1 (z)e x , (20)
where b and m are, respectively, the unknown constant temperature gradient and mass fraction gradient in the x-direction. For the stationary state, when the abovementioned assumptions are made and the corresponding boundary conditions are considered, we obtain the velocity, temperature and concentration fields: with

⎧ ⎪ ⎨ ⎪ ⎩ U (z) = -ψ 0 (2z -1), W 1 (z) = - ψ 0 Ra (2z -1), T 0 = -z + 1, C 0 = mx + g(z), (21)
⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ ψ 0 = 1 2 Raψm, g(z) = ψ 0 mLe(6z 2 -4z 3 -1) 12 + 1 -mA 2 -z, m = 5Leψ 0 Le 2 ψ 2 0 + 30 . ( 22
)
Then we have 2ψ 0 3 Le 2 + (60 -5RaLeψ)ψ 0 = 0 whose solutions are ψ 0 = 0 and ψ 0 = ± 1 2 √ 10RaLeψ-120 Le . The dimensionless and the dimensional mass fraction gradient are

m = ± √ 10RaLeψ -120 RaLeψ , m dim = √ 10q 1 q 2 HΔT -120 H 2 q 1 , ( 23
) where q 1 = Kgβ c /νD and q 2 = -C 0 (1 -C 0 )D T /D.
Figure 5 illustrates the evaluation of the dimensional mass fraction gradient m dim according to H and ΔT for a mixed water ethanol (60.88% by weight of water). It follows that for a fixed height of the cavity H, the mass fraction gradient will be so much greater than ΔT is great, while the H that maximizes the mass fraction gradient depends on the choice of ΔT , which is equal to H max = 16 q1q2ΔT .

Linear stability analysis of the unicellular flow

In order to analyse the stability of the unicellular flow, we introduce the perturbation of the vertical velocity component w, the perturbation of the vertical component of W , w 2 , the perturbation of temperature θ, the concentration c and the transformation η = cθ. The disturbances are developed in the form of normal modes. It is assumed that the perturbation quantities are sufficiently small so that the second-order terms may be neglected. The system of equations for the amplitudes can be written as

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ (Bσ + 1)(D 2 -k 2 )w + Rak 2 [(1 + ψ)θ + ψη] +Rv IkDw 2 b + ψ m +k 2 w 2 DT 0 + ψ DC 0 +Ik 3 w 2 [(1+ψ)θ + ψη] = 0, Ik(D 2 -k 2 )θ -σIkθ -wIkDT 0 + Uk 2 θ + bDw = 0, Ik Le (D 2 -k 2 )(η + θ) -σIk(η + θ) -wIkDC 0 +Uk 2 (η + θ) + mDw = 0, (D 2 -k 2 )w 2 + k 2 [(1 + ψ)θ + ψη] = 0. ( 24
)
The corresponding boundary conditions are given as

⎧ ⎪ ⎪ ⎨ ⎪ ⎪ ⎩ w = 0, θ = ∂η ∂z = 0, w 2 = 0 for z = 1, w = 0, ∂θ ∂z = ∂η ∂z = 0, w 2 = 0 for z = 0. ( 25
)
The linear stability equations were solved using the 4th-order Galerkin method:

⎧ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎨ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎪ ⎩ w = N n=1 a n (1 -z)z n , θ= b 1 z 2 -1 + b 2 z 3 -1 + N n=1 b n+2 z 2 -1 z n+1 , η = d 1 + d 2 z 2 - 2 3 z 3 + N n=1 d n+2 (1 -z) 2 z n+1 , w 2 = N n=1 e n (1 -z)z n .
(26) For the value of ψ and Le studied, the critical Rayleigh number leading to stationary bifurcation is always higher than the one leading to Hopf bifurcation. So, in this study we focus on the critical number related to the Hopf bifurcation. The results are illustrated in table 2 for the case ψ = 0.3, Le = 100, = 0.5, A = 10 and B = 10 -6 , the vertical vibrations have a stabilizing effect on convective motions and increase the value of the critical thermal Rayleigh number Ra c . Thus the vibrations can be used to maintain the unicellular flow and allow species separation over a wider range of Rayleigh numbers. It should be mentioned that vibrations reduce the critical wave numbers k c and the Hopf frequency ω c . This means that vibrations can also be used to decrease the number of convective cells at the transition from the unicellular flow to the multicellular flow.

2D Numerical simulations

The averaged equations (eq. ( 6)) with the associated boundary conditions were solved using the finite element method. The influence of vibrations on the onset of convection was investigated for a cell of aspect ratio A = 10 for Le = 100 and = 0.5 (B is fixed to 10 -6 ). It was observed that the critical parameters of the bifurcations differed very little between the case A = 10 and the case of a cell of infinite horizontal extention. A structured mesh 150 × 30 was used for the finite element method for A = 10. For the onset of stationary convection, the results for Le = 100 are presented in figs. shows the streamlines for Rv = 0. In this case, the flow is multicellular and we cannot use the horizontal cell to separate the components of a binary mixture. Figure 6(b) shows the streamlines for the same values of all the parameters but with vibrations. It can be observed that vibrations modify the structure of the flow from multicellular to unicellular. Similar results were obtained for ψ = 0.01 (fig. 7).

Conclusions

The two-dimensional thermosolutal Soret-driven convection under mechanical vibration was studied analytically and numerically. The vibration was in the limiting range of high frequency and small amplitude and its direction was taken parallel to the temperature gradient. The geometries considered were an infinite horizontal layer and a confined cavity. The linear stability analysis of equilibrium and unicellular solutions were studied for the positive value of ψ. The equilibrium solution was found to lose its stability via a stationary bifurcation. For the long-wave mode, the Rayleigh number was obtained as Ra = 12 Leψ . It was observed that the vertical vibrations have a stabilizing effect on convection. The vibrations could be used to decrease the value of the separation factor leading to the unicellular flow ψ uni , allowing separation in the binary mixture in a horizontal cell for a broad range of binary mixtures. The dimensional mass fraction gradient m dim was defined according to H and ΔT . It followed that, for a fixed height of the cavity H, the mass fraction gradient increases with ΔT , while the height that maximizes the mass fraction gradient is H max = 16 q1q2ΔT . The expression of the analytical relation is S = A • √ 10RaLeψ-120 RaLeψ giving that the separation value is identical to the one obtained by Elhajjar et al. [9] for a cell with horizontal walls maintained at constant but different temperatures.

Fig. 1 .

 1 Fig. 1. A rectangular cavity, H is the height and L is the width, the low horizontal plate is submitted to constant uniform heat flux and the top one is maintained to a constant temperature while no mass flux is imposed

Fig. 2 .Fig. 3 .

 23 Fig. 2. Effect of vertical vibration on the onset of convection for Le = 100, = 0.5 and B = 10 -6 .

1 1039Fig. 4 .

 14 Fig. 4. Effect of the vertical vibration on the onset of the unicellular separation ratio beyond which the flow at the onset of convection becomes unicellular for = 0.5 and B = 10 -6 .

Fig. 5 .

 5 Fig. 5. Curve of dimensional mass fraction gradient m dim for a mixed water ethanol (60.88% by weight of water), for (β c = 0.2, K = 10 -9 , C0 = 0.5, D = 4 • 10 -10 , DT = 10 -12 , ν = 3 • 10 -6 ).

  6 and 7. For ψ = 0.02, without vibration (Rv = 0), the critical parameters Ra c = 5.78, k c = 1.00 are obtained from the linear stability analysis. For the same value of ψ but with vibrations (Rv = 10), we obtain Ra c = 6, k c = 0 from the linear stability analysis, so the critical wave numbers is zero, which means that the flow at the onset of convection is unicellular. To confirm this result, we used the direct numerical simulation to study the case ψ = 0.02 for a value of Ra close to the critical value (Ra = 6) first without vibrations (Rv = 0) and then with vibrations (Rv = 10).

Fig. 6 .Fig. 7 .

 67 Fig. 6. Streamlines for Le = 100, ψ = 0.02, Ra = 6. (a) Rv = 0 (without vibration); (b) Rv = 10. A multicellular flow is obtained at the transition from the equilibrium solution for Rv = 0, whereas a unicellular flow is obtained for Rv = 10.

Figure 6 (

 6 Figure 6(a) shows the streamlines for Rv = 0. In this case, the flow is multicellular and we cannot use the horizontal cell to separate the components of a binary mixture. Figure6(b) shows the streamlines for the same values of all the parameters but with vibrations. It can be observed that vibrations modify the structure of the flow from multicellular to unicellular. Similar results were obtained for ψ = 0.01 (fig.7).

Table 1 .

 1 Comparison of the critical Rayleigh and the wave number for Le = 100, = 0.5 and B = 10 -6 .

	ψ	Rv = 0	Rv = 10	Rv = 50
		Rac	kc	Rac	kc	Rac	kc
	0	27.10 2.33 30.35 2.12 40.43 1.66
	0.001 23.56 2.20 26.55 1.97 35.29 1.48
	0.002 20.75 2.09 23.49 1.83 30.98 1.31
	0.004 16.59 1.90 18.88 1.57 24.16 0.99
	0.005 15.03 1.81 17.11 1.45 21.41 0.83

Table 2 .

 2 Effect of vibrations on the critical values of the Rayleigh number Ra c, wave number kc and frequency ωc associated with the transition from unicellular to multicellular flow for ψ = 0.3, Le = 100, = 0.5 and B = 10 -6 .

	Rv	kc	Rac	ωc
	0	5.25	16.87	1.55
	10	4.98	17.09	1.44
	30	4.54	17.39	1.27
	50	4.20	17.55	1.14
	70	3.94	17.59	1.04

This research is supported by the CNES (Centre National d' Études Spatiales).

Author contribution statement

AA and AM conceived of the presented idea and supervised the findings of this work. SO and AM developed the theory and performed the computations. AA and AM verified the analytical methods and encouraged SO to compare numerical simulations and analytical results. All authors discussed the results and contributed to the final manuscript.


R Céolin

Josep Ll Tamarit
email: josep.lluis.tamarit@upc.edu

Fullerene Solvate

J Ye

M Barrio

Ph Negrier

N Qureshi

I B Rietveld

J Ll

Tamarit

Orientational order in the stable buckminster fullerene solvate C60.2CBr2H2

niveau recherche, publiés ou non, émanant des établissements d'enseignement et de recherche français ou étrangers, des laboratoires publics ou privés.

Introduction

Fullerene C 60 is known to be moderately soluble in a large number of solvents, from small quasi-tetrahedral molecules to aromatic ones. Moreover, C 60 has a strong tendency to form co-crystals with solvent molecules, i.e. solvates .

In many cases the intermolecular interactions between C 60 and the guest solvent molecules are strong enough to form a co-crystal and weak enough to allow for the wellknown orientational disorder of the C 60 and solvent molecules. It has been demonstrated that due to the relatively weak intermolecular interactions, most solvates are not stable in air and lose solvent molecules through evaporation resulting in disorder within the crystalline structure of desolvated C 60 fullerene [26].

The orientational disorder in solvates can complicate the solution of structures by X-ray diffraction [27] and only in a few cases, C 60 molecules were found to exhibit orientational order as for example in the recent case of the low-temperature solvate consisting of C 60 and cubane (C 8 H 8 ), which orientational order is claimed to be caused by topological molecular recognition between the convex surface of C 60 and concave cubane [28,29]. Often, structural refinement is only possible by allowing for several rotational orientations of the C 60 molecule.

For example in certain cases, the solvate structure could be refined using site occupancy factors with two distinct orientations for C 60 , revealing an obstructed molecular rotation [30,31]. Moreover, it is often seen that the solvent molecules have different orientations within the C 60 interstices. It is clear therefore that the moderate interaction strength between the C 60 molecules and its guest molecules is highly tunable, which makes its solvates an interesting subject of study.

In this work, the solvate formation between C 60 and dibromomethane, CBr 2 H 2 , has been investigated. Dibromomethane is a relatively small, rigid molecule with C 2v symmetry and a significant dipole moment of 1.51 D. It is expected to cause anisotropic van der Waals interactions in line with other, previously investigated molecules [32][33][34]. We will demonstrate that in this case the intermolecular interactions are strong enough to form a cocrystal that is stable in air and in which the host and guest molecules are orientationally ordered.

Dibromomethane (CBr 2 H 2 , M= 173.835 g mol -1 ) was purchased from Aldrich with a purity higher than 99%. Fullerene C 60 (M = 720.64 g mol -1 ) was purchased from TermUSA (purity > 99.98%). Both compounds were used as received.

Solid C 60 was dissolved in CBr 2 H 2 at room temperature in screw-cap tubes that were subsequently stored for several months in the dark. Black, bright, polyhedral crystals formed at the surface of the colorless solution. The morphology of the crystals (Figure 1 was corrected following the recommended procedure by external calibration using the cubic phase of Na 2 Ca 3 Al 2 F 4 and cubic spline fitting [35]. Peak position determinations were carried out after pseudo-Voigt fitting. Crystals were gently crushed and introduced in a Lindemann capillary (0.5 mm diameter) rotating perpendicularly to the X-ray beam during the experiments to improve the average over the crystallite orientations.

A liquid nitrogen 700 series Cryostream Cooler from Oxford Cryosystems (U.K.) was used to control the temperature during the measurements. For the temperature dependent measurements the heating rate in between data collection was 1.33 K min -1 and the sample temperature was equilibrated for at least 10 min after which the X-ray diffraction profiles were recorded isothermally for at least 60 min. Longer acquisition time patterns, 9 h, were recorded for structural and refinement purposes at room temperature. Pawley fits, indexing, structure solution, and Rietveld refinement were performed with the Materials Studio Program (Version 5.5) from Accelrys (San Diego, CA) [START_REF]Materials Studio) version 5.5[END_REF]. An independent FullProf structural analysis was carried out to verify the refinement results.

Results

An X-ray diffraction pattern was obtained from the crystals (see Figure 1) extracted from the surface of the colorless solution at room temperature. After the crystals were left at room temperature in contact with air for several days, the same X-ray patterns were obtained.

It can therefore be assumed that the solvate is stable in air.

Systematic absences for the space group assessment were compatible with the centered monoclinic space group C2/m. To determine the structure of the solvate, atom coordinates of the individual molecules C 60 and CBr 2 H 2 were used from the literature and fitted in the obtained diffraction patterns. For the Rietveld refinement, a rigid-body constraint was used for CBr 2 H 2 . The position and orientation of the molecules were refined with a single overall isotropic displacement parameter and preferred orientation using the March-Dollase formula [START_REF] Dollase | [END_REF]. The refinement results are depicted in Figure 2, together with the experimental pattern and the difference between the refined pattern and the experimental one. The values of the final Rietveld refinement are summarized in Table 1.

Surprisingly, the C 60 molecules appeared to be orientationally ordered in the current solvate. For confirmation, an independent Rietveld refinement of the structural model was carried out with the FullProf program. In a first step, the C 60 molecule was modeled as a spherical shell with the scattering density of 60 carbon atoms homogeneously distributed on a sphere with a radius of 3.59 Å, which represents a freely rotating molecule. The CBr 2 H 2 molecule was modeled as a rigid body with the atom coordinates taken from Podsiadło et al. [32]. The local coordinate system was set with the 2-fold rotation axis of the CBr 2 H 2 molecule pointing along the local z axis (the local coordinate system is defined by x||a, y||b and z||(a x b) with θ=ϕ=χ=0), which therefore permits to rotate the molecule around its 2-fold rotation axis by the angle χ independently of its orientation in the unit cell. A first refinement clearly shows that the CBr 2 H 2 molecule is oriented with its 2-fold rotation axis parallel to the monoclinic b axis, i. e. θ=ϕ=90 degrees.

In a second step, the C 60 molecule was modeled as a rigid body with its 2-fold rotation axis along the local z axis. A clear improvement of the refinement is apparent between the spherical-shell model with the agreement factor, R F , equal to 8.0 and the rigid-body model with R F = 4.2, which strongly suggests that the C 60 molecule is not freely rotating. The resulting orientation is such that its 2-fold rotation axis is aligned with the monoclinic b axis just like the CBr 2 H 2 molecule (i.e. θ=ϕ=90 degrees).

Rotating the two molecules each around their individual 2-fold rotation axes corresponding to the angle χ for both molecules exhibits a strong influence on the agreement factor, in particular for dibromomethane, however also C 60 has a clear single minimum within one π periodicity (see Figure 3). The results confirm the orientational order of dibromomethane with the same alignment of its rotation axis in relation to the monoclinic cell. A dynamic rotation of the CBr 2 H 2 molecule around its 2-fold axis can be ruled out due to the very strong dependence of the angle χ on R F (see Figure 3), which demonstrates that the molecule must have a fixed orientation, i.e. with the Br ligands pointing towards the longer diagonal of the a-c plane (see Figure 4). The final structural parameters obtained from the last refinement were virtually the same as those obtained by means of the Materials Studio package. In the case of C 60 the results are not as clear-cut as for CBr 2 H 2 , but the difference between the global minimum and the local minima is convincing enough to expect C 60 to be mainly stuck in a single orientation, even if a stepwise arrested rotation cannot be excluded due to the many local minima. A few crystals together with a small quantity of mother liquor were subjected to thermogravimetric analyses. After the recorded mass reached a constant value at room temperature, the sample was assumed to be mother liquor free. It was subsequently heated to 550 K (Figure 5) with a constant heating rate of 2 K min -1 . Experimental mass loss was found to be ca. 32%, i.e. very close to the expected value of 33% ( = 2×173.8 / 1067.7) for the C 60 :CBr 2 H 2 = 1:2 molar ratio. The TG curve in Figure 5 clearly demonstrates that the desolvation process takes place in two overlapping steps. It could indicate the formation of a solvate with a lower stoichiometry or the existence of a desorption process following the desolvation of the 1:2 solvate.

The same procedure as for the TG measurements was followed for the DSC measurements (Figure 5). Again, it can be observed that two thermal effects are convoluted as the endothermic peak (T onset = 326 K) clearly consists of two overlapping peaks. In fact, the derivative of the mass loss with respect to the temperature (gray curve in Figure 5, which results from the TG measurements) mimics the DSC signal even if a small shift in temperature exists. Grey curve corresponds to the derivative of the weight loss with respect to the temperature.

The dashed red arrow indicates the mass loss at the minimum of the derivative (see text). The endothermic effect recorded by DSC and associated to the complete desolvation can be separated into peaks a and b, which respectively represent enthalpy changes of 73.9 J g -1 and

29.1 J g -1 related to the initial mass of the solvate.

To investigate the two convoluted peaks in the desolvation process more closely, X-ray diffraction measurements were carried out on crystals in a closed capillary as a function of the temperature. Figure 6 shows several patterns between room temperature and 385 K, which is near the end of the second endothermic process. Patterns at intermediate temperatures reveal that the Bragg peaks of the monoclinic solvate remain clearly present up to 330 K above which the diffraction pattern appears of desolvated C 60 partially amorphized by crystal defects and stacking faults as commonly observed after desolvation of C 60 solvates [16,21,22,39]. It indicates that the existence of a second solvate with a smaller stoichiometry is not likely.

Discussion

The desolvation process, which can be written as

C 60 •2CBr 2 H 2 (s) → C 60 (s) + 2 CBr 2 H 2 (g) (1)
, occurs in two steps as demonstrated by DSC as well as thermogravimetry by peaks a and b in Figure 5. The enthalpy change of the process is found to be 103 J g -1 of solvate (110.0 kJ mol - 1 of solvate) or 316 J g -1 of solvent (54.9 kJ mol -1 of solvent) by summing over peaks a and b as a whole. The total enthalpy is higher than the sublimation enthalpy of CBr The bimodal thermal effect associated with the desolvation process was not caused by the formation of a second solvate as demonstrated by the X-ray experiments as a function of the temperature in Figure 6, because impure fcc C60 is observed in the temperature range related to peak b and no additional peaks of a possible new solvate have been observed. Because the desolvation process must occur first, the second step will be vaporization of the released solvent. Thus, after the desolvation process, CBr 2 H 2 may first liquefy and remain temporarily adsorbed to the C 60 followed by evaporation. This can be summarized with the following two steps:

C 60 •2CBr 2 H 2 (s) → C 60 (FCC) + CBr 2 H 2 (l, adsorbed) (2)
C 60 (FCC) + 2 CBr 2 H 2 (l, adsorbed) → C 60 (FCC) + 2 CBr 2 H 2 (g)

Peak a, which has an enthalpy change of 73.9 J g -1 (78.95 kJ mol -1 ) of initial solvate, consists of the complete destruction of the solvate network, step (2), and the evaporation of about 1 mol of solvent, half the step (3). Thus, the enthalpy change related to step (2) is found by subtracting the vaporization enthalpy of one mole of solvent (i.e. step 3) from 78.95 kJ mol -1 of initial solvate. This leads to ΔH(2) = 32.8 kJ mol -1 of initial solvate. The enthalpy change related to peak b (29.1 J g -1 of initial solvate) belongs to 83.7 % of the remaining sample, which coincidentally has a 1:1 mol ratio (one mol of solvent per mol of C 60 ). The enthalpy can therefore be expressed as 34.77 J g -1 of remaining sample or as 31 kJ mol -1 of desorbing and evaporating solvent. This value is of the same order as the enthalpy needed to desorb similar organic solvents from graphitized thermal carbon [START_REF] Avgul | Chemistry and Physics of Carbon[END_REF].

The X-ray diffraction results and the Rietveld refinement reveal a monoclinic (space group C2/m) crystal structure with overall orientationally ordered C 60 and CBr 2 H 2 molecules. The refinement results appear to point to orientationally ordered molecules, both for C 60 and CBr 2 H 2 in the solvate structure. This is quite unusual, because in general in solvates with small halogen-methane or -ethane derivatives, the C60 molecules, and often the solvent molecules too, are found to exhibit orientational disorder. For example, solvent molecules possessing C 2v symmetry, such as CBr 2 Cl 2 or CBr 2 (CH 3 ) 2 , were found to be orientationally disordered in solvates with a hexagonal structure, whereas the packing of such solvates is very similar to the present one, as illustrated in Figure 7A.

The overall orientational order in the CBr 2 H 2 solvate cannot be explained by the dipole moment of the solvent of 1.51 D, because CBr 2 (CH 3 ) 2 possesses a similar value of 1.64 D in contrast to that of Br 2 CCl 2 of 0.2 D, which both exhibit orientational disorder in their respective solvates. To take a closer look at the host-guest interactions, excess volumes can be studied. Such volumes are defined as the difference between the measured volume of a solvate (defined as V unit cell /Z, i.e. 1701.5/2=850.75 Å 3 ) and the sum of the molecular volumes of the C 60 and solvent molecules from their respective pure structures. For C60, the molecular volume is 710 Å 3 from its FCC structure and the molecular volume of Br 2 CH 2 was determined as 95.175 Å 3 from the structure by Kawaguchi et al. [33]. Accordingly, the excess volume is found to be negative: 850.75 -(710 + 2×95.175) =-49.6 Å 3 , which is relatively high and demonstrates the strong interaction between the C 60 and solvent molecules.

The packing coefficient, η, of the present solvate can be determined as the van der Waals volumes of the solvate molecules divided by the measured solvate volume V unit cell /Z. The van der Waals volume of C 60 is known to be 526 Å 3 (the volume of a sphere with a ca. 5 Å radius)

and that of the CBr 2 H 2 molecule has been determined through the Kitaigorodsky method [START_REF] Kitaigorodskii | Mixed Crystals[END_REF] as 66.

Conclusions

The hexagonal solvate C 60 •2CBr 2 H 2 , which is stable in air, has been structurally and thermodynamically characterized. The C 60 molecules possess a hexagonal base structure, which is deformed and becomes monoclinic due to the presence of solvent molecules.

Notwithstanding the similarity of CBr 2 H 2 with other solvent molecules that form C 60 hexagonal solvates, in particular those with the same C 2v molecular symmetry as CBr 2 Cl 2 and CBr 2 (CH 3 ) 2 , the CBr 2 H 2 solvate exhibits overall orientational order for both the C 60 and the CBr 2 H 2 molecules, rarely seen in the other solvates. The orientational order is consistent with the solvate's stability, which has a high negative excess volume and a high desolvation enthalpy. All these physical properties demonstrate the strong host-guest interactions.

Although the correct space group symmetry of the solvate is monoclinic (C2/m), the overall packing is very similar to the hexagonal packing found in many other solvates.

It is clear that the present solvate with CBr 2 H 2 exhibits overall orientational order due to strong interactions between C 60 and its solvent molecule. This stands in clear contrast to the recent case of orientational order in the solvate mentioned in the introduction of C 60 and cubane, which was claimed to be caused by topological molecular recognition between the convex surface of C 60 and the concave cubane [28].

  ) was examined by means of scanning electron microscopy (JEOL-7100F) with a scanning voltage of 20 kV. Differential scanning calorimetry (DSC) and Thermogravimetry (TG) were performed under nitrogen flux at a 2K min -1 scanning rate using the Q100 and Q50 from TA Instruments (New Castle, DE, USA), respectively. Sample masses were weighed in at a precision of 0.01 mg.

Figure 1 .

 1 Figure 1. Scanning electron microscopy (SEM) photograph of C 60 •2CBr 2 H 2 crystals

Figure 2 .

 2 Figure 2. Experimental (red circles) and calculated (black line) diffraction patterns along with the difference profile (blue line) and Bragg reflections (vertical bars) of the monoclinic C2/m space-group of the C 60 •2CBr 2 H 2 solvate at room temperature. The inset corresponds to the scale for the data between 40 and 60° 2θ.

Figure 3 .I

 3 Figure 3. The crystallographic agreement factor, R F , as a function of the rotation angle around the 2-fold rotation axis of C 60 (blue dots) and of CBr 2 H 2 (red dots). A clear optimal

Figure 4 .

 4 Figure 4. The molecular positions in the crystal structure of C 60 •2CBr 2 H 2 approximately along the [100] direction (left panel) and along the [010] direction (right panel).

Figure 5 .

 5 Figure 5. Differential scanning calorimetry curve (black, left axis) and thermogravimetry curve (red points, right axis) as a function of the temperature for the C 60 •2CBr 2 H 2 solvate.

Figure 6 .

 6 Figure 6. X-ray patterns of the C 60 •2CBr 2 H 2 solvate as a function of the temperature. The pattern of pure FCC C 60 is shown at the top for reference purposes.

  The structure is characterized by alternating planes consisting of either C 60 or guest molecules stacked along the c axis. The CBr 2 H 2 molecule is located at [0, 0.2710(2), 0.5] in the unit cell, while the position [0, 1/3, 1/2] corresponds to the prismatic void [1/3, 2/3, 1/2] of the hexagonal parent structure of C 60 with space group P6/mmm. In fact, the monoclinic b axis in the C-centered cell is approximately 2 × cos(30 ○ ) times the hexagonal a axis. Therefore, due to the identical packing and void filling, the structure of C 60 •2CBr 2 H 2 can be regarded as a distorted hexagonal one.

Figure 7 .

 7 Figure 7. (A) Measured, weighted unit-cell volumes (V/Z) of C 60 solvates against the van der Waals volume of its solvent molecule and (B) the packing coefficient as a function of the van der Waals volume of its solvent molecule. In (B), the values for the monoclinic (higher value) and orthorhombic (lower value) of 2C 60 •3CS 2 solvate have been added. Dashed lines are guides to the eyes.

Table 1 . Crystal structure information and Materials Studio results of the Rietveld refinement for the C 60 •2CBr 2 H 2 solvate

 1

	Chemical Formula	C 60 •2CBr 2 H 2
	M / g•mol -1	1068.31
	2θ-Angular Range	7 -80 o
	Space group	C2/m
	a /Å	9.9001±0.0004
	b / Å	17.446±0.001
	c / Å	10.1013±0.0004
	α / o	90
	β / o	102.769± 0.002
	γ / o	90
	V/Z / Å 3	1701.5±0.4
	Z (Z')	2(1)
	Temperature	293 K
	D x / g•cm -3	2.085± 0.001
	Wavelength (Cu Kα1)	λ=1.5406 Å
	2θ-shift (zero correction)	0.0264 ± 0.0012
	Profile Parameters
	Na	0.509 ± 0.015
	Reliability Parameters

  2 H 2 of 46.2 kJ mol -1 (obtained by adding the melting enthalpy, 9.2 kJ mol -1 , and the vaporization enthalpy, 37 kJ mol -1[40]). It implies that CBr 2 H 2 and C 60 attract each other relatively strongly in the C 60 •2CBr 2 H 2 solvate. The TG-curve indicates that the first peak is accompanied with a weight loss of 16.3 % of the initial mass. This corresponds to the loss of about one of the two moles of CBr 2 H 2 from the solvate C 60 •2CBr 2 H 2 as 0.163 × 1068.31 g mol -1 = 174.1 g mol-1 , which is very close to the molar mass of the solvent of 173.8 g mol -1 .

  3 Å 3 . With these values, a packing coefficient has been found of 0.773, which is higher than the packing coefficient of close packed C 60 (η = 0.74) and that of the monoclinic (C2/m) low-temperature ordered phase of CBr 2 H 2 at 183 K (η = 0.70). It demonstrates that the interaction between C 60 and CBr 2 H 2 is strong confirming the likelihood of an overall

		1000	(A)		CCl 2 Br 2 CCl 2 (CH 3 ) 2			0.78	CS 2	CBr 2 H 2	(B)
	) V/Z(Å 3	900	CCl 3 H CCl 3 (CH 3 ) CBr 2 ClH CCl 3 Br CBr 3 H CCl 4	η	0.74 0.76	CCl 2 Br 2 CBr 2 ClH CCl 3 H CBr 3 H CCl 2 (CH 3 ) 2
			CBr 2 H 2						CCl 4	CCl 3 Br
			50	70	90	110	0.72	50	70	90 CCl 3 (CH 3 ) 110
			VvdW-Solvent (Å	3 )			VvdW-Solvent (Å	3 )
	orientational order of the solvate molecules. Figure 7B illustrates the high packing coefficient
	for the CBr 2 H 2 solvate compared to other hexagonal solvates containing C 60 and similar
	solvent molecules. The only comparably sized solvate that exhibits orientational order for C 60
	and the solvent molecule too, the low-temperature orthorhombic solvate with CS 2 has a
	similar packing coefficient, ca. 0.78, as the present solvate (cf. Figure 7B). The corresponding
	high-temperature monoclinic solvate of CS 2 on the other hand, which exhibits orientational
	disorder, possesses a considerably lower packing coefficient, ca. 0.76 (see Figure 7B) [29].

Acknowledgements

This work has been supported in part by the Spanish Ministry of Science and Innovation (grants FIS2014-54734-P) and the Catalan Government (grant 2014 SGR-581).