J.-C Géminard

D Bouraya

H Gayvallet

Thermal conductivity associated to the bead-bead contact decorated by a liquid bridge An experimental study based on the response of a chain subjected to thermal cycles

Keywords: . 66.70.+f Nonelectronic thermal conduction and heat-pulse propagation 05.Rm Porous materials, granular materials -91.60.Ki Physical properties of rocks and minerals, thermal properties

We report an experimental study of the heat transport along a chain of macroscopic beads: at one end, one bead is periodically heated and we record the resulting temperature variations of another bead as a function of position and time. The experimental results show that the chain behaves like a high-order low-pass filter. The measurement of the associated cut-off frequency makes possible to determine accurately the resistance of the bead-bead contact: we make use of the experimental setup for studying the effects of liquid bridges decorating the bead-bead contacts, which provides important clues for understanding the thermal properties of partially wet granular matter.

Introduction

The transport of electricity, sound or heat in a dense assembly of solid particles addresses interesting and puzzling questions. Even if granular materials are widely used in industry, one is still far from understanding entirely how the transport coefficients relate with the physics of the contacts between the particles and the geometry of the contact network (texture). For instance, for over a century [1,2], electrical transport in metallic powders has generated interest. These systems exhibit fascinating properties, such as highly non-linear and hysteretic electrical conductivities as well as extreme sensitivity to electromagnetic waves [3]. The propagation of soundwaves in granular materials also exhibits very interesting features; nonlinearities such as a slight dependence of the transmission amplitude on the frequency of the source at low vibration amplitude and an extreme sensitivity to the position of each grain [4]. The propagation of soundwaves is so sensitive to the texture that it is likely to be altered by the vibration itself or by small changes in the temperature. However, the propagation of non-linear waves in a onedimensional chain of beads under static stress has been experimentally proven to follow the prediction of Hertz theory as long as plastic deformations of the grains at the contacts are not involved [5].

In the same way, the transport of heat in granular materials addresses puzzling questions which are relevant to be answered as thermal properties of granular matter are of great practical importance. Indeed, the most common materials used to insure thermal as well as acoustic insulation in buildings consist in assemblies of solid particles or fibers. For instance, glasswools are efficient thermal insulators that also present the advantage to be low-cost, imputrescible and fire resistant as they are produced mainly from recycled glass and sand [6]. Moreover, because of their ductility, they can absorb any unevenness of the substrate. Another important application worth to be mentionned here is the use of granular materials as buffers in nuclear-waste deposits [7]. Experimental and theoretical studies of the heat transport in granular materials have already been reported [2,[START_REF] Kling | [END_REF][10][11]. They pointed out the crutial role played by the physics of the contacts between the grains. Accordingly, humidity has been experimentally proven to increase significantly the thermal conductivity as water tends to decorate the contact points [12].

Whereas this latter study considered the thermal properties of the bulk granular material, involving thus the texture of the material, there does not exist, to our knowledge, any experimental study dedicated to the thermal properties of the decorated contacts alone. In the present article, we report an experimental study of the thermal conductivity of a one-dimensional chain of beads in the same spirit as preceding studies of electric conductivity [3] and sound propagation [5] that made it possible to analyze the physics of the contacts in absence of texture effects. The response of the system to temperature cycles shows that a periodic medium can be an efficient thermalinsulator in the sense that, for instance, the chain is very efficient in filtering temperature changes as it acts as a high-order low-pass filter. We determine the corresponding cut-off frequency in order to obtain reliable measurements of the thermal resistance associated to the beadbead contact. We make use of the experimental setup for studying the effects of liquid bridges decorating the beadbead contacts.

The content of the manuscript is organized as follows: in the section 2, we introduce the principle of the experiment and describe the experimental setup and procedure. In the section 3, we report our experimental results obtained in the cases of dry or wet contacts and compare to theoretical predictions. Finally, we conclude in the section 4.

Experimental setup and procedure

The experiment consists in heating periodically one end of a linear chain of beads and in recording the resulting temperature of the n th bead as a function of time t. In the next section 2.1, we describe the experimental setup whereas the section 2.2 is dedicated to the procedure used to analyze the thermal response of the system to the heat injection (Fig. 1).

Experimental setup

The chain consists of 10 centimetric ball-bearing steelbeads (AISI 304, diameter d = 1 cm, from Marteau & Lemarié) aligned with the help of three PTFE plates [see the chain cross-section in the inset (Fig. 1)] that insure small thermal contact with the remaining part of the experimental setup which consists of a stainless steel frame. The contact between the beads is insured by pushing (static axial force F ) the chain at one end against a fixed steel cylinder located at the other end. A force sensor (Entran, ELA-B2E-10KN) is used to measure F to within 10 N (F < 10 kN).

In order to inject heat at one end of the chain, the first bead is equiped with a small heating wire (Constantan, diameter 100 µm, typical resistance r ≃ 1 Ω) placed at its center. We impose the periodic voltage U (t) = U 0 1 + sin 2πνt from a home-made linear power-amplifier driven by a function generator (Stanford Research Systems, DS345). A first multimeter (Keithley, 2001) is used to measure the voltage difference U (t) across the heat source (Fig. 2). In our experimental conditions, the frequency ν ranges from 1/3600 to 1/30 Hz whereas the maximum voltage 2U 0 is of about 3 V. The power P (t) injected in the first bead thus writes:

P (t) = P 0 3 4 + sin 2πνt + 1 4 sin 4πνt + 3π 2 (1)
where P 0 = 2U 2 0 /r (we neglect here the variation of r with the temperature). We point out that the bead is continuously heated with the mean power 3 4 P 0 and that the variation with time of the heating power P (t) contains two The chain consists in 10 centimetric steel-beads align with the help of three PTFE plates (inset). In order to insure the contact between the beads, a static force F is applied at one end (bead 10). The first bead (bead 1) is heated by means of a small resistive wire located at its center: the heating voltage is measured by means of a first multimeter. In order to obtain the local temperature of the chain, the bead n (n = 6 in the sketched situation) is equiped with a Pt-sensor which resistance is measured with the help of a second multimeter. In addition, the temperature of the surrounding air, between beads 7 and 8, is measured by means of a second Pt-sensor (not sketched in this figure).

harmonics. We shall make use of this latter property for studying simultaneously the response of the system at the three frequencies 0, ν and 2ν. U Fig. 2. Experimental T (2) and U vs. time t. We report, for ν = 1/3600 Hz, the temperature T (2) of bead 2 and the voltage U as functions of time t during 2 hours in the permanent regime (the total duration of the experiment is 10 hours, corresponding to 10 periods of the excitation). The mean temperature is about 2.5 K above the temperature Ta of the surrounding air (δT (2) 0 ≃ 2.5 K). The amplitude of the temperature variations at ν is about 2.21 K whereas it is about 0.46 K at 2ν (dry contacts, F = 100 N).

The bead n (the index n refers to the position of the bead in the chain) is equiped with a temperature sensor (Pt100, from Heraeus) located at its center. A second multimeter (Keithley,196), used in 4-wires configuration, measures the resistance R s of the sensor so that the local temperature of bead n, T (n) (t), is known to within about 10 -3 K. As we do not regulate the overall temperature of the experimental setup, we measure, in addition, the outside air temperature, T a , by means of a second temperature sensor (Pt100) located in air and connected to a third multimeter (Keithley,196).

The function generator and the three multimeters are connected to a computer through a IEEE interface. For given experimental conditions (position n of the sensing bead, applied force F , dry or wet contacts), a simple Croutine automatizes the experiment: for a chosen set of heating powers P 0 and frequencies ν, the voltage U (t) and the temperature T (n) (t) are measured (sampling frequency 1 Hz) and saved to the hard-drive (Fig. 2). Analysis of the data is performed afterwards (Igor Pro 4, Wave-Metrics, Inc.).

Transfer function

From equation (1), if the system has a linear response, we expect the temperature T (n) (t) to exhibit the same spectral components as P (t), this is to say, harmonics at frequencies 0, ν and 2ν. From the raw data (Fig. 2), we extract five relevant quantities; the amplitudes, δT between the bead n and the outside atmosphere (0 frequency). To do so, we interpolate the experimental temperature-difference

∆T (n) (t) ≡ T (n) (t) -T a with ∆T (n) (t) = δT (n) 0 +δT (n) 1 sin (2πνt + φ (n) 1 ) +δT (n) 2 sin (4πνt + φ (n) 2 ).
(

) 2
The origin of time t in equation ( 2) is accurately obtained by interpolating U (t) with a sine-function and by choosing the origin such as U (0) = U 0 and dU dt (0) > 0. We first check that the response of the system is linear: the amplitude of each of the harmonics, δT

(n) 1(,2) , scales like U 2
0 (thus, is proportional to the heating power) whereas the corresponding phase, φ (n) 1(,2) , does not depend on U 0 (Fig. 3).

Thus, from the quantities δT

(n) 0 , δT (n) 1 , δT (n) 2 , φ (n) 1 and φ (n)
2 , we can determine experimentally the power-to-(temperature of the bead n) transfer-function, n) , as follows: a first set of data points is obtained from the first harmonics at the frequency ν. From the first harmonics, we get

H (n) (ν) ≡ H (n) e jθ (
H (n) = δT (n) 1 and θ (n) = φ (n) 1 .
In equation ( 1), we notice that the power injected in the first harmonics at ν is 4 times larger than that injected in the second harmonics at 2ν, and that there is phaseshift 3π 2 between them. Thus, from the same experiment, the second harmonics provides another set of data points as

H (n) = 4δT (n) 2 and θ (n) = φ (n) 2 -3π
2 at the frequency 2ν. Note that, by construction, H (n) has the dimension of temperature and is proportional to P 0 . Strictly speaking, we should divide H (n) by P 0 to get the power-to-(temperature of the bead n) transfer-function, but, as we vs. U 2 0 . As an example, we report the amplitude, δT

( 10 -3 K ) 0.4 0.3 0.2 0.1 0.0 U 0 2 ( V 2 ) 3.0 2.0 1.0 0.0 φ 1 (3) ( rad ) δT 1 (3) φ 1 (3)
1 , and phase, φ

1 , measured for bead 3 as functions of the square heating-voltage U 2 0 . The amplitude of the harmonics δT

(n) 1(,2) is proportional to U 2 0 and the phase Φ (n) 1(,2)
does not depend on U0, which proves that the response of the system is linear.

will not change P 0 in what follows, we will report H (n) as the temperature defined above. In figure 4, we report the typical experimental result for H (2) (ν). We observe that the data evaluated from the first and second harmonics superpose down to temperature variations of amplitude about 10 -3 K.

θ (2) (rad) φ 1 (2) φ 2 (2) -3π/2 δT 1 (2) 4 δT 2 (2)
Fig. 4. Experimental H (2) and θ (2) vs. frequency ν.

As an example, we report the amplitude, H (2) , and the phase, θ (2) , as functions of the frequency ν. The data evaluated from the first or the second harmonics are marked with different symbols. We observe that the responses estimated from the first and second harmonics respectively superpose down to temperature variations of amplitude about 10 -3 K (the lines are only guides for the eye).

We shall explain in the section devoted to the experimental results (section 3) how the response of the chain can be extracted from these data, provided we are able to model the response of the experimental setup.

Experimental procedure

The experimental procedure is as follows: for chosen experimental conditions (applied force F , dry or wet contacts), in order to obtain the response of the system, we first place the bead containing the temperature sensor at a chosen position n. For a series of excitation frequencies ν, we record the temperature variation T (n) as a function of time t in the permanent regime. The response at a different position n is obtained after opening the experimental setup and changing the position of the sensing bead. The experimental procedure makes possible to determine H (n) (ν) for n ∈ [2,7] in the frequency range 1/ν ∈ [30, 3600] s. Measurements at frequencies smaller than 1/3600 Hz are, in principle, possible but note that we record the temperature variations for 10 periods, which then would lead to an experimental time larger than 10 hours. For frequencies larger than 1/30 Hz, the temperature variations become small (amplitude typically less than 10 -3 K) and the response of the system is determined with poor accuracy. We always make sure that the measurements are reliable by checking that the phase of the measured harmonics is locked (i.e. that the phase is not arbitrary with respect to that of the heating voltage) and exclude the data points that do not satisfy this condition.

Experimental results

In the following, we make use of the experimental setup for studying the thermal properties of the chain in different experimental conditions. In the next section 3.1, we consider the case of dry contacts in air: the analysis of the whole set of experimental data makes possible the modeling of the thermal response of the whole experimental setup. Then, section 3.2 is dedicated to the study of wet contacts obtained by addition of liquid droplets in the bead-bead contact-regions.

Dry contacts

First experimental observations

In figures 5 and 6, we present the whole set of experimental data obtained for the chain of steel beads with dry contacts and F = 100 N (the static applied force is small and only insures that the beads are at contact). At first sight, the chain behaves like a low-pass filter which order depends on the position n of the bead within the chain (the larger n is, the larger is the order of the filter). Nevertheless, we cannot immediately extract the caracteristic frequency of the filter without exploring in details the response of the experimental setup.

For instance, consider only the response of the system at zero frequency. Due to thermal losses, one would expect the mean temperature δT (accordingly H (n) (0)) to vanish (i.e. T (n) (t) → T a ) for n → ∞. However, even if we observe that H (n) (0) decreases almost exponentially when n is increased, the mean temperature-difference H (n) (0) surprizingly does not vanish for large n (Fig. 7). The exponential decay of H (n) (0) to a finite temperature-difference is the most obvious experimental manifestation of the nontrivial response of the experimental setup. It turns out that, the response at zero frequency is not enough for us to measure the thermal losses to the outside atmosphere: the whole set of experimental data is necessary for obtaining a correct model of the experimental situation, which we discuss in the next section. We observe that the temperature difference at zero frequency between the bead n and the atmosphere decreases almost exponentially when n increases but does not vanish far away from the heat source (large n).

Dynamical response

From the experimental data presented in figure 5, we can extract the dependance of the amplitude H (n) on the po-sition n at a given frequency ν (Fig. 8). We observe again that, even at non-zero frequency, the amplitude H (n) decreases exponentially when n increases but does not vanish for large values of n. We write:

H (n) (ν) = H (0) (ν) exp [-k ′ (ν)n] + H (∞) (ν) (3)
and report in figure 9 the experimental values of H (∞) (ν) [Eq. ( 4)] and k ′ (ν) = ℜ(k) [Eq. ( 5)] which can be accounted by the modeling of the experimental situation as follows (Fig. 10): We consider the beads as heat capacitors (heat capacity C) connected one to its two neighbors through thermal resistors (resistance R), and we assume that the thermal losses to the frame can be associated to the thermal resistance, β along the chain (n ≥ 2) and γ at bead 1. Detailed description of the model is proposed in the appendix. We observe that, for a given frequency ν, the amplitude H (n) does not vanish far away from the heat source (large n). The result holds true in the frequency range explored experimentally and we measure that the limit H (∞) decreases when the frequency ν is increased (Fig. 9).

H (∞) (ν) (K) k' H (∞)
Fig. 9. Asymptotic amplitude H (∞) and k ′ vs. ν. We observe experimentally that the asymptotic amplitude H (∞) decreases, according to equation ( 4), when the frequency ν is increased. In addition, k ′ = ℜ(k) obeys equation (5), which makes possible to caracterize the thermal conduction along the chain.

We observe a good agreement of the experimental data (Fig. 9) with

H (∞) (ν) ∼ 1 1 + 2πν ωc 4 (4)
for ω c = (0.8 ± 0.01) 10 -3 rad s -1 where ω -1 c ≡ C √ βγ from the model [see appendix, equation ( 18)]. In addition, we expect k ′ to be the real part of k, defined by

cosh (k) = 1 + 1 ω 0 τ c + j 2πν ω 0 . ( 5
)
where ω -1 0 ≡ RC/2 and τ c ≡ βC. From equation ( 5) and the experimental data, we get ω 0 τ c = 1.87 ± 0.03 and ω 0 = (7.6 ± 0.1) 10 -3 rad s -1 . Note that the first quantity, ω 0 τ c ≡ 2β/R, is the ratio of the resistance β associated to the energy losses to the resistance R of the bead-bead contacts whereas ω 0 ≡ 2/(RC) characterizes the transport of heat along the chain alone.

Quantitative analysis

In this section, we shall compare the experimental value of the contact resistance, which we can extract from the experimental ω 0 , with a theoretical modeling due to G. K. Batchelor and R. W. O'Brien [2].

First, one can find in the litterature the heat capacity c s = 0.5 J g -1 K -1 and the density ρ s = 8 g cm -3 of steel (AISI 304), which lead to C = 2.1 J K -1 for d = 1 cm. Thus, from the experimental value of ω 0 , one obtains the experimental value, R exp = (126 ± 1) K W -1 , of the resistance of the bead-bead contact.

Second, we can consider theoretically the heat flux Ψ ≡ (T (n) -T (n+1) )/R from bead n to bead (n + 1). According to the Hertz theory, due to the axial applied force F , the radius ρ of the contact region is given by [START_REF] Landau | Theory of elasticity[END_REF]:

ρ = F 1/3 3(1 -σ 2 )d 8E 1/3 (6)
where E = 2 10 11 Pa is the Young modulus and σ = 0.29 the Poisson ratio of the material the beads are made of. In our experimental conditions (F = 100 N), we find ρ = 0.12 mm. In addition, the heat conductivity of steel, k s = 16.2 W m -1 K -1 , is large compared the thermal conductivity of air, k a = 0.025 W m -1 K -1 . Denoting ζ = k s /k a and ξ = 2ζρ/d (in our experimental conditions, ζ ≃ 650 and, from the radius ρ of the contact region, ξ ≃ 15.6), we can write

2Ψ πk a d(T (n) -T (n+1) ) = H 0 (ζ) + H c (ξ) + ∆H m (ξ). ( 7
)
In equation ( 7), H 0 (ζ) stands for the dimensionless heatflux across the air layer in the case of a point-like contact.

Provided that ζ ≫ 1,

H 0 (ζ) = 2 log e k s k a + A -3.9 ( 8
)
where A is a constant of the order of unity which depends on the outer field conditions and is independant of ζ. The second contribution H c (ξ) in equation ( 7) accounts for the dimensionless heat-flux across the contact circle of finite radius ρ. In the limit ξ ≫ 1, H c (ξ) ∼ 2ξ/π. We point out that, from [2], H c (ξ) obeys this asymptotic behavior typically for ξ > 10. Finally, the last contribution to the heat flux in equation ( 7), ∆H m (ξ) accounts for the difference between the real flux across the air layer and the flux expected in the case of a point-like contact. In the limit ξ ≫ 1, ∆H m (ξ) ∼ -2 log e ξ but we note that this asymptotic behavior is expected for very large values of ξ. From the experimental value of ξ, we expect ∆H m (ξ) ≃ -0.5 [2]. Finally, from equation (7), we obtain the theoretical value R theo = 2546/(18.5 + A) K W -1 which agrees quantitatively with the experimental value R exp ≃ 126 K W -1 for A ≃ 1.7, which is of the order of unity. Thus, we conclude from the excellent agreement (in spite of the uncertainty in the value of A) between R exp and R theo that measurements of the frequency ω 0 provides reliable values of the contact resistance R.

In the next section 3.2, we make use of the experimental setup for analyzing the effect of liquid bridges decorating the regions of contact between the beads.

Wet contacts

As already mentionned in the introduction, humidity has been experimentally proven to increase significantly the thermal conductivity of an assembly of glass beads as water decorates the contact points [12]. Our experimental setup is particularly suitable for studying the effects of liquid decorating the contact points: indeed, a change of the physical conditions in the contact region is expected not to affect any of the model parameters (τ c , τ b , α as defined in the appendix) excepted the cut-off angular frequency ω 0 , which relates to the resistance of the contact. Thus, one can save experimental time and study the evolution of ω 0 alone by analyzing the response of the system at a given position n only.

Experimental procedure

We shall report the evolution of the thermal resistance R of the bead-bead contact as a function of the size of the liquid bridges that decorate the contact regions. As an example, we choose to study the response of the system at the bead 4: we know from the study of the dry system that the amplitude of the temperature variations is large enough for giving significant results and that the order of the 'filter' is also large enough for providing accurate measurements of the cut-off frequency ω 0 . This choice is anyway in a large manner arbitrary. We measure the diameter Φ of the liquid bridge on photographs of the bead-bead contact-region (Φ = 2.17 mm, the lateral size of the photograph is 1 cm).

Nonvolatile liquid bridges (identical along the chain) are created by placing a given amount of silicon oil (Rhodorsil, 47V5000) in each of the bead-bead contact-region. We expect, from the large viscosity of the liquid (ν = 5000 cP), the temperature gradient not to induce any significant convective flow and the transport of heat to be only due to diffusion. We image the system from side (Fig. 11), which makes possible to measure the diameter, Φ, of the liquid bridge (the variability along the chain remains less than 2%). We measure φ before and after measurement of the thermal response of the system in order to make sure that the bridges where equilibrated during the experimental time.

Results

We report in figure 12, the amplitude H (4) as a function of the frequency ν for the dry system and three different diameters Φ. The experimental results exhibit an excellent agreement with the theoretical expression (18) in which we use the experimental parameters, τ b = 318 s, τ c = 246 s, α = 2289 and P τ b /C = 13.3 K determined previously for the dry system and adjust only the value ω 0 of the cut-off frequency. From the experimental value of ω 0 , provided that the thermal inertia of the bridges remains negligible compared to that of the beads, we can easily determine ) Φ = 2.17 mm Φ = 1.75 mm Φ = 1.00 mm Φ = 0.00 mm Fig. 12. Amplitude H (4) vs. 2πν. The amplitude H (4) is affected by the presence of liquid bridges in the contact regions between the grains: we observe that, at a given frequency ν, H (4) increases with the diameter Φ of the bridges. The presence of liquid leads to a decrease in the thermal resistance R of the contact. The experimental data (symbols) are successfully accounted by the theoretical description of the system [lines, from equation ( 18) with different values of ω0]. We observe that the thermal resistance of the contacts decreases when the diameter of the liquid bridge is increased. The dashed line is a guide for the eye while the horizontal black line corresponds to the value, R∞ = 56.6 K W -1 , estimated for beads totally immersed in oil [from equations (7) and [START_REF] Landau | Theory of elasticity[END_REF] with A = 1.7]. The grey line corresponds to equation [START_REF] Kling | [END_REF]. the thermal resistance of the contact R which is found to decrease when the bridge diameter φ is increased (Fig. 13).

The theoretical description of the thermal resistance associated with a bead-bead contact of finite size decorated by a liquid bridge is, to our knowledge, not yet available. The extension of the model briefly summarized in section 3.1.3 to this latter case appears to be difficult and we will limit the discussion of the experimental results to a few comments.

First, we can obtain an estimate of the thermal resistance R ∞ , one would measure for beads immersed in silicon oil, and compare this latter value to the minimum of R reported in figure 13. The thermal conductivity of the silicon oil k so = 0.16 W m -1 K -1 . In this case, ζ ≃ 101 and ξ ≃ 2.4 (for F = 100 N). Replacing k a by k so in equation ( 8), we get H 0 (ζ) ≃ 5.33 + A, which accounts for the heat transport through the fluid layer around the contact (A is again a constant of the order of unity). Moreover, we know from reference [2] that, for a value of ξ as small as 2.4, the value 2ξ/π ≃ 1.5 largely overestimates the real value of H c (ξ), which is associated to the conduction through the contact circle of radius ρ, and that the correction ∆H m (ξ) remains much smaller than unity. Thus, we can neglect H c (ξ) and ∆H m (ξ) in equation (7), which leads to R ∞ ≃ 398/(5.33 + A) K W -1 . If we assume that the value A = 1.7 remains unchanged when the beads are immersed in air or in silicon oil, we get R ∞ ≃ 56.6 K W -1 (Fig. 13). We observe that the resistance of the contact R approaches R ∞ when the diameter Φ of the bridge is increased. Even if R remains about 30% larger than R ∞ for Φ/d ≃ 0.3, we consider that the agreement between our experimental data and the theoretical estimate of R ∞ again reinforce the validity of our measurements. Moreover, we are tought by the model that, in the case of beads immersed in oil, the heat is mainly tranfered through the liquid layer around the contact region and not through the circle of contact (indeed, H c (ξ) ≪ H 0 (ζ)). We expect the conclusion to hold true in the case of a large liquid bridge decorating the contact region. This povides us with an important clue for understanding the smooth variation of R with the diameter φ of the bridge.

We could expect from equations ( 7) and ( 8) and from the ratio air-to oil-conductivity, R to suddenly decrease when the liquid is introduced in the contact region. However, we observe experimentally that R decreases almost linearly when Φ is increased. This feature may be qualitatively understood as follows: in the case of dry beads, about half of the heat is transfered through the real beadbead contact, the second half passing through the air gap around the contact in spite of the small radius ρ of the contact circle. Think now about what would append, if the radius ρ was slightly increased by δρ. We would mainly expect H c (ξ) to increase linearly with δρ as H c (ξ) ∝ ρ, ∆H m (ξ) to remain negligible, and H 0 (ζ) not to change. Note that to increase ρ by δρ is thus equivalent to decorate the contact with a steel bridge of diameter Φ = 2(ρ + δρ). The corresponding variation δR of the resistance would then be given by δR = -Rδρ/ρ {from equation (7), δ[log e (R)] = -δ[log e (H c )]}. Thus, introduction in the contact region a small quantity of oil, which conductivity is smaller than that of steel, is expected to only continuously decrease the resistance of the contact and not to produce any discontinuous change in R. If we assume that the addition of oil produces the same effect as the addition of steel excepted that the transport of heat occurs in oil instead of in steel, we can roughly estimate

δR R = - k so k s Φ -2ρ 2ρ (9)
which leads to dR dΦ ≃ -5.2 K W -1 mm -1 . As the presence of the liquid bridge modifies also the temperature field in the bead, leading to a spreading of the temperaturefield lines and, thus, to reduced thermal resistance. The slope given by equation ( 9) is expected to overestimate the experimental value for small Φ, which is observed experimentally (Fig. 13). Obviously, due to assumptions which remain open to criticism, this quantitave estimate must be considered with caution. However, we consider that the reasoning presented above is qualitatively correct for explaining the smooth dependancy of the resistance R on the diameter Φ, when the contact region is decorated by a liquid bridge.

Conclusions

First, we would like to point out several flaws of our experimental situation; due to long thermal times (especially 2π/ω 0 ≃ 15 mn) the experimental times are long [for instance, more than 10 hours are needed for measuring H (n) (ν) at five relevant frequencies ν]; the study of R as a function of the applied force F or of the nature of the surrounding gas, which is in principle possible, would require a huge experimental effort. Indeed, as at least the times τ b and τ c would depend on F or on the nature of the surrounding gas, characterization of the whole experimental setup would be necessary. The analysis could not be limited to the response of the system at a given position n in the chain, which would increase significantly the experimental time. Note however that the bevavior of H (n) (ν) as a function of n at one given frequency ν would be, in principle, enough for determining ω 0 . Such studies remain to be done although we already attempted, before drawing our conclusions, to measure the resistance of the contacts in air, nitrogen, carbon dioxide and helium and obtained sensible qualitative-results. However, these studies were limited to the thermal response at a given position n as a function of the frequency ν, which unfortunately turned out not to be sufficient to determine R afterwards.

In conclusion, we have shown that the study of the heat transport along a chain of beads is a pertinent way to access the thermal properties of the contact. The experimental method is applied successfully to the case of dry contacts in the sense that the measured resistance R can be accounted by a classical model by Batchelor et al [2]. We made use of the experimental setup for studying the effect of liquid bridges decorating the contact regions and observed a smooth decrease in the thermal resistance when the diameter of the bridge is increased. Even if a theoretical description of this latter experimental situation is not yet available, simple arguments, which extend the model of the dry contact, allow to understand, at least qualitatively, the experimental results. This work provides important clues for understanding the transport of heat in partially wet granular materials. n increases and tends to a finite value for large values of n. These experimental observations are accounted by the following modeling of the experimental system (Fig. 10).

The thermal conductivity of the chain can be accounted as follows : we consider each bead as a thermal mass, C, in contact with its two neighbors (excepted bead 1) through identical thermal resistors, R. We assume that the thermal losses through the PTFE plates, which are used to align the beads, can be described as the effect of a single thermal resistor, β, which connects locally each bead to the remaining part of the experimental setup. Let us denote T (n) (t) the temperature of the bead n at time t. We can thus write, taking ω 0 ≡ 2/(RC) and τ c ≡ βC,

∂T (n) ∂t = ω 0 2 T (n+1) + T (n-1) -2T (n) - 1 τ c T (n) -T f (x n ) (n > 1) (10)
where T f (x n ) denotes the temperature of the experimental setup at the position x n ≡ (n -1)d of the bead n (diameter d). We will assume that the chain is long enough to be considered as semi-infinite: the boundary condition far away from the heat source then only writes T (n+1) = T (n) , (n → ∞). We will further assume that the contact between the first bead with the metallic part, that constrains the chain along its axis, can be described by a simple termal resistor, γ. The first bead is heated by injecting the heating power P (t) so that the energy balance leads to, taking τ b ≡ γC, ∂T

Let us now assume that the experimental setup, apart from the PTFE plates and the chain, behaves like a semiinfinite rod (diffusion coefficient D) in contact with the outside atmosphere (at temperature T a ) through the conductance per unit length 1/η. Furthermore, we assume that the thermal contact between the beads and the rod along the chain is bad enough for the temperature of the rod not to be affected by the corresponding heat sources (i.e. β is large and the diffusion in the rod is fast). As a consequence, we can write :

∂T f ∂t = D∆T f - 1 τ l (T f -T a ). ( 12
)
where τ l = ηC f is the characteristic time associated with the thermal losses to the atmosphere (C f is here the heat capacity per unit length of the metallic frame). We will consider that the temperature profile T f (x, t) only results from energy injection at x = 0 from bead 1 (we thus assume that γ ≪ β). Writing the energy balance at the contact between bead 1 and the rod leads to :

T (1) -T f (0) γ = -DS∇T f (0) ( 13
)
where S stands for the cross section of the rod. After having written, the general set of equations ( 10) and ( 12) and the boundary conditions (11) and (13), we consider the harmonic response of the experimental setup as a function of the frequency ν = ω/(2π) by writing :

T f (x, t) = T 0 exp(jωt -k f x), (14)
T (n) (t) = T (n) 0 exp(jωt), (15)
and P (t) = P (ω) exp(jωt). ( 16)

From equation ( 12), we get k 2 f = η+jω D whereas the boundary condition (13) leads to T and T 0 are obtained as functions of P by using the second boundary condition (11).

In the limit |k f d| ≪ 1, which is well satisfied in our experimental conditions, equation ( 17 We point out that the diffusion coefficient in the rod, D, only appears, in this final expression proposed for |T (n) 0 |, in the coefficient α. We are in the limit |nk f d| ≪ 1, which leads to |k f d| ≪ 1 when we consider the limit n ≫ 1 (the diffusion in the rod is fast) but, even in this limit, α ≡ γDSk f can remain finite. Note that, from equation (18), the amplitude T (n) 0 of the temperature oscillations decreases exponentially when n is increased and reaches, for n large, the value, T (∞) 0 which magnitude remains finite as observed experimentally [Fig. 9, Eq. ( 4)].

Fig. 1 .

 1 Fig. 1. Sketch of the experimental principle.The chain consists in 10 centimetric steel-beads align with the help of three PTFE plates (inset). In order to insure the contact between the beads, a static force F is applied at one end (bead 10). The first bead (bead 1) is heated by means of a small resistive wire located at its center: the heating voltage is measured by means of a first multimeter. In order to obtain the local temperature of the chain, the bead n (n = 6 in the sketched situation) is equiped with a Pt-sensor which resistance is measured with the help of a second multimeter. In addition, the temperature of the surrounding air, between beads 7 and 8, is measured by means of a second Pt-sensor (not sketched in this figure).

2 ,

 2 of the two harmonics at the frequencies ν and 2ν and the mean temperature difference δT (n) 0

Fig. 3

 3 Fig. 3. Amplitude δT(3) 1 and phase φ(3) 1

Fig. 5 .Fig. 6 .

 56 Fig.5. Amplitude H (n) vs. frequency ν.The amplitude H (n) decreases faster with increasing frequency ν for larger n (F = 100 N, dry contacts, the lines are only guides for the eye).

Fig. 7 .

 7 Fig.7. Experimental H (n) (0) vs. bead position n. We observe that the temperature difference at zero frequency between the bead n and the atmosphere decreases almost exponentially when n increases but does not vanish far away from the heat source (large n).

  (ν) = (0.09 ± 0.01) K k' = 1.04 ± 0.02 ν = 1/3600 Hz

Fig. 8 .

 8 Fig.8. Amplitude H (n) vs. bead position n. We observe that, for a given frequency ν, the amplitude H (n) does not vanish far away from the heat source (large n). The result holds true in the frequency range explored experimentally and we measure that the limit H (∞) decreases when the frequency ν is increased (Fig.9).

Fig. 10 .

 10 Fig. 10. Sketch of the experimental situation.

Fig. 11 .

 11 Fig. 11. Photograph of a liquid bridge.We measure the diameter Φ of the liquid bridge on photographs of the bead-bead contact-region (Φ = 2.17 mm, the lateral size of the photograph is 1 cm).

Fig. 13 .

 13 Fig. 13. Resistance R vs. bridge diameter Φ.We observe that the thermal resistance of the contacts decreases when the diameter of the liquid bridge is increased. The dashed line is a guide for the eye while the horizontal black line corresponds to the value, R∞ = 56.6 K W -1 , estimated for beads totally immersed in oil [from equations(7) and (8) with A = 1.7]. The grey line corresponds to equation[START_REF] Kling | [END_REF].

=

  T 0 (1 + α) with α ≡ γDSk f . In addition, equation(10) imposes : (n -1)k f d] 1 + ω 0 τ c [1 -cosh (k f d)] + jωτ cwith cosh (k) = 1 + 1 ω0τc + j ω ω0 . The constants T (0) 0

ω 0 τ b [ 1 -

 1 [α + j(1 + α)ωτ c ] exp [-(n -1)k] / (1 +jωτ c ) [α + j(1 + α)ωτ b ] exp (-k)][α + j(1 + α)ωτ c ]

Acknowledgments

The authors would like to thank C. Tassius and E. Freyssingeas for careful reading and useful criticism of the manuscript.

Appendix

Modeling of the experimental situation

We observe experimentally that the temperature oscillation at a given frequency ν decreases exponentially when


P.-E Roche

Applicability of Boussinesq approximation in a turbulent fluid with constant properties

Keywords: numbers: 47.27.Te Turbulent convective heat transfer, 44.25.+f Natural convection, 47.55.P-Buoyancy-driven flows; convection Boussinesq approximation, Natural convection, Rayleigh-Bénard, Turbulence

The equations of motion describing buoyant fluids are often simplified using a set of approximations proposed by J. Boussinesq one century ago. To resume, they consist in assuming constant fluid properties, incompressibility and conservation of calories during heat transport. Assuming fulfilment of the first requirement (constant fluid properties), we derive a set of 4 criteria for assessing the validity of the two other requirements in turbulent Rayleigh-Bénard convection. The first criterion α∆ ≪ 1 simply results from the incompressibility condition in the thermal boundary layer (α and ∆ are the thermal expansion coefficient and the temperature difference driving the flow). The 3 other criteria are proportional or quadratic with the density stratification or, equivalently with the temperature difference resulting from the adiabatic gradient across the cell ∆ h . Numerical evaluations with air, water and cryogenic helium show that most laboratory experiments are free from such Boussinesq violation as long as the first criterion is fulfilled. In ultra high Rayleigh numbers (Ra > 10 16 ) experiments in He, one of the stratification criteria, scaling with α∆ h , could be violated. This criterion garanties that pressure fluctuations have a negligible influence both on the density variation and on the heat transfer equation through compression/expansion cycles. Extrapolation to higher Ra suggests that strong violation of Boussinesq approximation could occur in atmospheric convection.

I. INTRODUCTION

The buoyant motion of fluids combines the non-linear complexity of hydrodynamics with additional degrees of freedom brought in by the thermodynamics. The exact equations for the buoyant flows are therefore very difficult to tackle directly. In 1903, Joseph Boussinesq introduced a set of approximations and derived simplified equations for buoyant motion [START_REF] Boussinesq | Théorie analytique de la chaleur[END_REF]. His three key approximations are : 1) the fluid's properties are constant, 2) the flow is incompressible (except for a linear temperature dependence of density yielding the buoyant term), 3) the heat equation is not coupled to the flow mechanical energy ("conservation of calories") Today, Boussinesq approximations are widely used to describe a large number of flows and to perform numerical simulations. A recent review presents an overview of works dedicated to the applicability of Boussinesq approximation and to its mathematical justification [START_REF] Zeytounian | [END_REF]. Oberbeck's name is occasionally associated with Boussinesq's to refer to the set of approximations presented above. As a side remark, it is worth noting that his convection model preserves compressibility and neglects all the time derivatives with respect to space derivatives [3].The resulting differential equations and physics significantly differs from the ones obtained by Boussinesq's.

To be able to predict if Boussinesq's approximations are applicable, specific quantitative criteria are necessary for each type of flow. The present work derives such criteria for Rayleigh-Bénard (RB) turbulent convection [4,5], assuming constant properties of the fluid (first approximation satisfied). A vast literature is dedicated to this first approximation, in particular for high Prandtl numbers in connection with mantle convection (for example see [6,7,8,9] and references within). For the intermediate Prandtl numbers, recent progresses have been reported in the turbulent regime (for example see [10,11,12,13,14] and references within).

Boussinesq second approximation ("incompressibility") and third one ("conservation of calories") consist in neglecting nearly half of the terms present in the exact equations of motion. When the magnitude of all the terms can be estimated using typical scales of the flow, the applicability of the approximation can be predicted. In laminar RB convection, temperature and velocity gradients smoothly spread across the cell, which eases the choice of typical scales in the flow. A set of criteria for the applicability of Boussinesq approximation have been proposed in this low Rayleigh number (Ra) limit (for example see [15,[START_REF] Tritton | Physical Fluid Dynamics[END_REF] and ref. within). In turbulent RB convection (high Ra), choosing the correct typical scales is less straightforward due to the flow complexity. In the present work, we attempt to do so using various experimental results and simple modelling. We derive a new set of 4 criteria for the validity of Boussinesq approximation in a turbulent RB cell. In our derivation, we prioritise high Rayleigh numbers Ra and intermediate Prandtl numbers (P r) fluids (say P r ∼ 0.5 -30), as this range of parameters is relevant to many geophysical flows of air and water, industrial flows and high Ra experiments with air, water, SF 6 and helium. Using these criteria, we investigate the validity of Boussinesq approximations in experimental set-ups with water, air and helium. In particular, we show that one criteria is much more stringent than the others in most experiments, at least below Ra ∼ 10 16 typically. At higher Ra, a second criteria could become relevant. We also examine if deviations from Boussinesq approximation can explain the apparent incompatibility of heat transfer measurements among very high Ra experiments. We find that it is not the case.

The next section recalls the exact equations of motion and singles out the different terms which have to be neglected to obtain Boussinesq's equations. In the third section, the magnitudes of these terms are compared with the retained terms by modelling turbulent convection. This yields a set of 4 criteria for the validity of Boussinesq approximation, which are summarised and discussed in the fourth section. The last section is a short conclusion.

II. THE BOUSSINESQ EQUATIONS

The temperature, pressure and density fields are decomposed into a constant (subscript "0"), a stratification (subscript "h") and a residual (θ, p and δρ) contribution :

T = T 0 + T h + θ (1a) P = P 0 + P h + p (1b) ρ = ρ 0 + ρ h + δρ (1c)
By definition, T 0 is the average of the top and bottom plates' temperatures. T h is the linearised adiabatic gradient term defined by Eq.(2a). The pressure P 0 + P h and density ρ 0 + ρ h fields satisfy the hydrostatic balance Eq.(2b) and thermodynamic equilibrium Eq.(2c) of a static stratified fluid with a temperature distribution T 0 + T h (z). At cell's mid-height (z = 0), we set ρ h (0) = 0 and P h (0) = 0.

T h = - αgT 0 c p z (2a) dP h dz = -g (ρ 0 + ρ h ) (2b) 1 ρ 0 + ρ h dρ h dz = -α dT h dz + χ dP h dz = α 2 gT 0 c p -χg (ρ 0 + ρ h ) (2c)
where α, g, c p and χ are the isobaric thermal expansion coefficient, the gravity acceleration, the isobaric heat capacity per unit mass and the isothermal compressibility coefficient.

The residual terms θ, p and δρ account for fluctuations and for the inhomogeneities which are not captured by the stratification terms T h , P h and ρ h . For example, fluctuations are generated by turbulence and inhomogeneities result from the temperature forcing and confinement by the cell side walls. The potential temperature difference ∆ and the adiabatic-gradient temperature difference ∆ h across a convection cell of height h are defined as :

∆ = θ(-h/2) -θ(h/2) ∆ h = T h (h/2) -T h (-h/2) = αgT 0 h c p
(the notation ∇ 2 will be used for the Laplacian operator to avoid confusion with the temperature difference ∆). The mass, momentum and energy conservation equations are :

∇v = - 1 ρ Dρ Dt (4a) ρ Dv Dt = -∇P + ρg + F (4b) ρc p DT Dt = k∇ 2 T + ρΦ + αT DP Dt (4c)
where D/Dt is the particle derivative and where F and Φ account for the viscous force and corresponding dissipation and where k is the thermal conductivity.

In this section, we assemble -step by step-a set of self consistent criteria which are sufficient to derive the Boussinesq equations from the conservation equations above. Since the successive approximations to derive Boussinesq equations are intricated, it is convenient to assume from the beginning that the density field is nearly homogeneous, or more precisely that the following two criteria are fulfilled (we omit the absolute values for better readability) :

ρ h ≪ ρ 0 (5a) δρ ≪ ρ 0 (5b)
A. A preliminary thermodynamics relation

Using the generalised Mayer relation and defining γ = c p /c v , we note that :

1 -γ -1 = α 2 T 0 ρc p χ = α∆ h ρ 0 gχh (6)
This equation relates the density stratification resulting from two different contributions : α∆ h is the thermal expansion associated with the adiabatic-gradient and ρ 0 gχh is the compressibility associated with the hydrostatic pressure. Since γ is of order one or higher (critical point vicinity), this relation shows that the two contributions are comparable in strength.

B. Momentum conservation

Using the definitions Eq.( 1b) and (1c), the momentum conservation equation Eq.(4b) can be written :

ρ Dv Dt = -∇p + δρg + F -∇P h + (ρ 0 + ρ h )g
By definition of P h , the last two terms cancel out. Anticipating on the solenoidal condition ∇v = 0, the friction is restricted to the shear term η∇ 2 v which appears in Navier-Stokes equation (zero bulk viscosity). Using the criteria Eq.( 5), Boussinesq's momentum equation :

ρ 0 Dv Dt = -∇p -ρ 0 αθg + η∇ 2 v (7)
is recovered if the extra criterion is fulfilled :

δρ ≃ -ρ 0 αθ (8)
In the above equation, δρ is approximated at first order. We are assuming implicitly that higher order terms don't contain qualitative features which are significant for the flow dynamics. Such a "significant" feature is present -for example-in the viscous term (dissipation) and, consequently, it cannot be neglected with respect to the inertial term Dv/Dt, although the later is known to be larger in turbulent regimes.

C. Mass conservation

Using the definition Eq.(1c), the mass conservation Eq.(4a) can be written :

∇v = - 1 ρ Dρ Dt = - 1 ρ Dδρ Dt - v z ρ dρ h dz
Using the definition Eq.(2c), the criteria Eq.( 5) and (8), and the thermodynamics relation Eq.( 6), we get :

∇v ≃ α Dθ Dt + v z ρ 0 gχ γ
Boussinesq's approximation requires that the velocity field is solenoidal (∇v = 0), which is satisfied if each term in the right-hand side (RHS) of the previous equation can be neglected with respect to the typical velocity derivatives ∂v i /∂x i of the left-hand side (LHS) ∇v, where i ∈ {x, y, z} :

v z ρ 0 gχ γ ≪ ∂v i ∂x i (9a) α Dθ Dt ≪ ∂v i ∂x i (9b)

D. Energy conservation

Using the definitions Eq.( 1) and ( 2), the energy conservation equation Eq.(4c) can be written :

ρc p Dθ Dt -k∇ 2 θ = ρΦ + αT Dp Dt + v z gα [T 0 ρ -T (ρ 0 + ρ h )]
(10) Defining κ = k/ρc p as the thermal diffusivity, Boussinesq's thermal equation (a "conservation of calories" equation) :

Dθ Dt -κ∇ 2 θ = 0 ( 11
)
is recovered if each term of the RHS of the equation Eq.( 10) can be neglected. A sufficient condition for this is found using the criteria Eq.( 5) and ( 8) and writing that T 0 provides a reasonable order of magnitude of T . We find the new criteria :

Φ/c p ≪ Θ ′ (12a) αT 0 ρ 0 c p Dp Dt ≪ Θ ′ (12b) v z gα c p [T h + θ (1 + αT 0 )] ≪ Θ ′ (12c)
where

Θ ′ ∈ { ∂θ ∂t , v i ∂θ ∂xi , κ ∂ 2 θ ∂x 2 i } III.

CRITERIA OF APPLICABILITY OF THE BOUSSINESQ'S APPROXIMATION

In this section, we evaluate the criteria ( 5), ( 8), ( 9) and ( 12) for the special case of high Ra (turbulent convection) and intermediate P r.

To check the first criterion Eq.(5a), the equations (2c) and (2b) could easily be solved exactly but an order-ofmagnitude solution for ρ h is enough for our purpose. Using also Eq.( 6), we easily find that the criterion Eq.(5a) becomes :

ρ 0 ghχ γ ≪ 1 or α∆ h γ -1 ≪ 1 (13)
To estimate the other criteria, we discriminate between two regions of the cell where the flow phenomenologies are known to differ : the bulk and the boundary layers near the top and bottom plates. Each criterion should be fulfilled in the two regions. Before proceeding with these two regions, it is useful to present order-of-magnitude fits of various quantities in the cell.

A. Few orders of magnitude

The scalings of various flow quantities have been reported in the literature of turbulent convection, in particular in the Rayleigh-Bénard cells. In this subsection, our aim is to present a few simple fits that will be used to estimate and compare the orders of magnitude of various contributions. We are not suggesting here that the chosen power law fits have an underlying physical meaning. We neglect the inhomogeneities within the boundary layers and within the bulk of the flow, which result from finite size effects e.g. confinement of the large scale circulation.

Below Ra = 10 12 and for 0.7 < P r < 30, all reported N u(Ra) data are in reasonable agreement (for studies of both Ra and P r in this range, see for example [START_REF] Verzicco | [END_REF]18,19,20,21]). We will use the fit :

(Ra < 10 12 ) N u ≃ 0.065Ra 1/3 (14)
The 1/3 exponent chosen in this fit is slightly larger than the exponents fitted in most studies (which are closer to 0.30-0.32), even when the temperature drop in plates is compensated [22]. But this simple fit still provides good numerical estimates for our purpose. Above Ra = 10 12 , the reported N u(Ra) dependence fall into two sets : data from 2 experiments roughly follows the scaling of Eq.( 14) ( [23,24]) while the others exhibit a steep increase of the Ra exponent from 0.31 up to typically 0.38 [25,26,27,28].The enhanced-heat-transfer law can be fitted by : (10 12 < Ra < 10 14 ; P r ≃ 10) N u ≃ 0.015Ra 0.38 The apparent incompatibility between the two sets of data is not understood and motivates present research (see also the numerical results [29,30]). Although all the data showing a heat transfer enhancement (including [31]) are consistent with the "ultimate regime" predicted by R. Kraichnan [32], an extensive characterisation of this new regime is still missing. Therefore, we will not attempt to derive specific Boussinesq criteria for it in this paper. In the Discussion section, we will come back on this point to discuss very high Ra experiments.

A large scale circulation or "wind" is known to occur in Rayleigh-Bénard cells. The strength of the l.s.c velocity v lsc can be characterised by the fit :

Re = hv lsc ν ≃ 0.15Ra 1/2 P r -3/4 (15)
which is a compromise between studies from various groups [28,33,34,35,36,37,38] conducted with various operational definitions of Re for 0.7 < P r < 30, based on local velocity or turn over time measurements [39]. An alternative estimate of the typical velocity in the bulk of the flow can be obtained from the rms velocity fluctuations. Comparing the fit of Re (above) with rms velocity data from the literature [35,37,40,41] we note that the scaling of the root-mean-squared (rms) velocity is approximately the same as the scaling of the large scale circulation velocity, and that the magnitudes are comparable (rms velocity is typically 40% of v lsc ). We will therefore consider only a single characteristic velocity at large scales. We should mention that, in the literature, the Ra dependence of Re has been fitted with effective exponents as low as 0.40 [35] (and even 0.36 ± 0.05 at some location in a square cell [42] ) and as high as 0.65 (in the horizontal plane [37]).

Many consistent studies report scalings of the rms temperature fluctuations in the bulk of the cell versus Ra (for example see [23,24,43,44,45]) but, to our knowledge, only Daya and Ecke [40] reports both a Ra and P r dependence. The fit below is derived from Daya and Ecke original data. We should mention that the reported P r dependence is poorly fitted by a power law but that such a fit still provides a fair order of magnitude for our purpose. θ ⋆ rms is defined as the rms temperature fluctuations normalised by the temperature difference ∆ across the cell.

θ ⋆ rms ≃ 0.3Ra -0.13 P r -0.5 (16)
Finally, we will need an estimate of the time derivative ∂θ/∂t in the bulk of the flow. The references [46,47] report corresponding measurements over a large range of Ra (5.10 8 < Ra < 10 15 ) for Prandtl numbers of order 1-10 (see also [45]). They interpret a transition occurring near Ra ∼ 10 11 with a change in power laws. Nevertheless, within a factor 2, we can still re-fit all their data with the following power law, which is enough for our purpose (the Prandtl number dependence is unknown) :

∂θ/∂t ≃ 0.7Ra 0.4 ∆κ h 2
In the bulk of the flow, we expect advection to be the most effective heat transport mechanism, which implies ∂θ/∂t ≃ v∇θ. The relation shows that we can also estimate the time derivative from v∇θ, taking v ≃ νRe/h and using recent temperature-gradient measurements for 10 9 < Ra < 10 10 and P r ≃ 5.4 [48]. Averaging the measured gradient from the various locations inside the cell [48], we obtain the fit :

∂θ/∂t ≃ v∇θ ∼ v|∇θ| ≃ 70Ra -0.165 Re ∆ν h 2
It is satisfactory that both estimations are in reasonable agreement (magnitude and Ra exponent). We will favour the first fit because its has been validated over a wider range of Ra. It is convenient to define a dimensionless derivative θ ⋆ ∂t as :

θ ⋆ ∂t = ∂θ ∂t h 2 ∆κ ≃ 0.7Ra 0.4 x (17)
where x = (P r/5) X accounts for the unknown Prandtl dependence through the unknown parameter X. For intermediate P r, we will assume that x is of order 1.

B. Boundary layer region

The present knowledge of scaling in the boundary layers is not as extended as in the bulk, in particular due to the complex interplay between the thermal and the velocity boundary layers (for recent progresses, see [49,50,51]). Nevertheless, it is well established that each boundary layer concentrate nearly half of the total temperature drop ∆ over a typical distance of order λ = h/2N u. An internal velocity scale v bl can be obtained writing that total heat flux experiences a crossover within the thermal boundary layer from a molecular diffusion transport (right on the plate) to a convective transport at a distance λ from the plate, which gives v bl = κ/λ = 2N u.κ/h. An external velocity scale is also imposed on the thermal boundary layers by the large scale circulation present in the bulk of the cell v lsc = Re.ν/h. Using the estimates of N u and Re given previously, we note that v lsc > v bl . In the following, we will choose the velocity scale which gives the more stringent criteria ("sufficient conditions"). Using these temperature, length and velocity scales, we evaluate below the criterion (5b), ( 8), ( 9) and (12) .

Local thermodynamic equilibrium implies the differential equation :

dρ = -ραdT + ρχdP
Subtracting the stratification contribution (2c) from this local equilibrium, the criteria Eq.(5b) and ( 8) are found to be equivalent to :

αθ ≪ 1 (18a) χp ≪ αθ (18b)
In the boundary layer, we have θ ∼ ∆/2. The order of magnitude of p is assessed taking the divergence of Eq.( 7) and using the solenoidal property of the velocity field. We find a Poisson equation for the pressure :

∇ 2 p = -ρ 0 ∇ [(v∇) v] + αρ 0 g ∂θ ∂z (19a)
The first term of the RHS gives the usual dynamical pressure of order :

p ∼ ρ 0 2 (v bl 2 + v lsc 2 ) ∼ ρ 0 2 v lsc 2
and the second RHS term, a buoyant source of pressure variation, can be estimated by integration as :

p ∼ gρ 0 α∆λ 4
Substituting the expressions v lsc = νRe/h and λ = h/2N u and using the fits for N u and Re, we find that the buoyant contribution to the pressure is smaller than the kinematic contribution. Neglecting the former and using Eq.( 6), both criteria Eq.( 18) becomes, up to a factor 2 :

α∆ ≪ 1 (20a) P rRe 2 Ra α∆ h 1 -γ -1 ≪ 1 (20b)
The second criterion Eq.( 20b) is comparable or less stringent than Eq.( 13), so we can drop it.

In the boundary layer, the next criterion Eq.(9a) is estimated as :

vρ 0 gχ γ ≪ v λ
After simplification by the typical velocity v and using the expression of λ, we find that the above condition is less stringent than the criterion Eq.( 13). We will therefore ignore it too.

To evaluate the criterion Eq.(9b) in the boundary layer, the particle derivative of the temperature is estimated with a diffusion process (heat is indeed transferred between the plates and the fluid by diffusion in the boundary layers) :

Dθ Dt ≃ ∆/2 λ 2 /κ (21)
and the space derivative of the velocity is estimated from the typical internal velocity as v bl /λ. The external velocity v lsc is not chosen because it would give a more loose condition. Substituting the expressions of Dθ/Dt, v bl and λ in the criterion Eq.(9b) gives criterion Eq.(20a) again.

We now consider the three criteria Eq.( 12) associated with the energy equation. The three RHS term Θ ′ of these inequalities can be estimated with the characteristic scales used above. All turn out to have the same typical magnitude, as expected from the chosen boundary layer model :

Θ ′ ≃ ∂θ ∂t ≃ v i ∂θ ∂x i ≃ κ ∂ 2 θ ∂x 2 i ≃ κ∆ 2λ 2
We estimate the LHS of the inequalities Eq.( 12) using v lsc as the typical velocity ("sufficient condition").

Φ c p ν c p v lsc λ 2 αT 0 ρ 0 c p Dp Dt < αT 0 ρ 0 c p ρ 0 v lsc 2 /2 λ/v lsc v z gα c p [T h + θ (1 + αT 0 )] v lsc gα c p ∆ h 2 + ∆ 2 (1 + αT 0 )
Substituting the expressions of v lsc (Re) and λ(N u) in the previous equations, we find up to prefactors of order unity :

∆ h T 0 × 2P r 2 Re 2 Ra ≪ 1 (22a) α∆ h × P r 2 Re 3 2RaN u ≪ 1 (22b) ∆ h T 0 ∆ h ∆ + αT 0 + 1 × P rRe 4N u 2 ≪ 1 (22c)
The second condition (22b) is comparable or more stringent than Eq.( 13) and will replace it. The αT 0 term in the parenthesis of Eq.(22c) gives a condition comparable or less severe than (22b) and can therefore be dropped. The unity term in the parenthesis gives a condition that can be compared to Eq.(22a) and it turns out to be typically up to 10 times more stringent for low Ra (Ra ∼ 10 7 ), comparable for Ra ∼ 10 13 , and less stringent above. We will come back on these 2 conditions later. The ∆ h /∆ term of Eq.(22c) gives the condition :

∆ 2 h T 0 ∆ × P rRe 4N u 2 ≪ 1 (23)

C. Bulk region

In the bulk, for the same reasons as in the boundary layers, the criteria Eq.(5b) and ( 8) are equivalent to :

αθ ≪ 1 χp ≪ αθ
The first condition is weaker than the criteria Eq.(20a) since the magnitude of temperature fluctuations in the bulk is significantly smaller than the total temperature drop across the cell (θ ⋆ rms ≪ 1). To assess the second condition, the magnitude of the pressure p is estimated using the Poisson equation Eq. (19). The buoyant contribution is estimated by analogy with electrostatic, αρ 0 g∂θ/∂z being the pressure source with a negligible remote contribution from the boundary layers[63] and a local contribution from the bulk.

p ρ 0 2 v 2 + αρ 0 g| ∂θ ∂z |l 2
where l is a correlation length of the temperature gradient ∂θ/∂z in the bulk. The buoyant contribution | ∂θ ∂z |l 2 can be estimated from the temperature rms fluctuation ∆θ ⋆ rms , the typical temperature derivative ∂θ/∂z i , and using the inequality

| ∂θ ∂z |l 2 < | ∂θ ∂z |hl ∼ | ∂θ ∂z |h ∆θ ⋆ rms |∂θ/∂z| ∼ h∆θ ⋆ rms
Estimating the typical velocity as v = νRe/h and using the fit of Re and θ ⋆ rms , we find that the buoyant contribution to the pressure variation is comparable or smaller than the kinematic term for Ra > 10 8 . We will ignore it. A sufficient condition for the criteria χp ≪ αθ in the bulk is therefore :

α∆ h 1 -γ -1 × P rRe 2 2Raθ ⋆ rms ≪ 1 ( 24
)
This new criterion is comparable to Eq.(22b). We choose the new one because it is derived from tighter estimates.

To estimate the criteria Eq.(9a), we take v z ∼ νRe/h as the typical vertical velocity. The typical derivative ∂v i /∂x i is estimated assuming that the velocity fluctuations in the bulk of the flow are driven by the large scale circulation on a large scale h and follows a Kolmogorovlike cascade. Then, the energy injected into cascade is related to the viscous dissipation Φ, which is proportional to the square of velocity spatial derivatives ( [START_REF] Monin | Statistical Fluid Mechanics: Mechanics of Turbulence[END_REF][START_REF] Grossmann | [END_REF]) :

v 3 lsc h = Re 3 × ν 3 h 4 ≃ Φ ∼ ν ∂v i ∂x j 2 (25)
from which we get the estimate :

∂v i ∂x i ∼ Re 3/2 ν h 2
Using those estimates and Eq.( 6), it is straightforward to show that Eq.(9a) in the bulk is less stringent than Eq.( 13). We can therefore ignore this condition. As a comment, the LHS of Eq.(9a) has been estimated using the magnitude of the velocity at large scales while the RHS has been estimated at small scales where the derivative is known to be larger. The solenoidal condition Eq.(9a) is nevertheless still accounted for at large scale only, as can be checked estimating the RHS as v lsc /h (one gets Eq.( 13)).

To estimate the next criterion Eq.(9b), we need an upper bound for the particle derivative Dθ/Dt. We can consider the relaxation by molecular diffusion of a λ-thick plume of excess temperature ∆/2. The result is the same as in the boundary layer : Eq.( 21). The criterion Eq.(9b) becomes in the bulk : 2N u 2 P rRe 3/2 α∆ ≪ 1 Using typical fits Re(Ra), N u(Ra), we find that the previous inequality is less severe than (20a) and can be ignored.

We now consider the three criteria Eq.( 12) associated with the energy equation. Boussinesq thermal equation Eq. (11) states that the advection of temperature inhomogeneties is only affected by the molecular thermal diffusion. In the boundary layers, we showed that these two transport processes have the same efficiency. In the bulk, we expect advection to become the dominant process, except in some localised regions (on the thermals interface) where diffusion could be as efficient as in the boundary layers. To estimate the three inequalities Eq.( 12), we will require that each LHS term is significantly smaller than the advection term v∇θ ≃ ∂θ/∂t = θ ⋆ ∂t ∆κ/h 2 of the RHS. It should be noted that the molecular diffusion term in the thermal equation provides a cut-off (or "dissipation") mechanism for temperature fluctuations at small scales and that we are not aiming to compare the magnitude of this cut-off term with the LHS terms of Eq.( 12).

To estimate the LHS of the first thermal criterion Eq.(12a), Eq.( 25) provides an order of magnitude of the viscous dissipation Φ which is found consistent with the Direct Numerical Simulation (DNS) results of Verzicco [41]. The criterion Eq.(12a) in the bulk becomes Eq.(26a).

To estimate the next criterion Eq.(12b), we assume that pressure and temperature fluctuates with similar typical time scales. The derivatives Dp/Dt and ∂θ/∂t of Eq.( 12b) are replaced with their integrated value p ∼ ρ 0 v 2 /2 and θ ∼ ∆θ ⋆ rms . We find Eq.(26b). To estimate the LHS of the third thermal criterion Eq.(12c), we take v z ∼ Re.ν/h, θ ∼ ∆θ ⋆ rms and T h ∼ ∆ h /2 to find Eq.(26c).

∆ h T 0 × P r 2 Re 3 Raθ ⋆ ∂t ≪ 1 (26a) α∆ h × P rRe 2 2Raθ ⋆ rms ≪ 1 (26b) ∆ h T 0 ∆ h 2∆θ ⋆ rms + αT 0 + 1 × P rReθ ⋆ rms θ ⋆ ∂t ≪ 1 (26c)
Using the typical orders of magnitude of Re, ∆θ ⋆ rms and θ ⋆ ∂t , we now compare these three thermal conditions Eq.( 26) with previous ones. The first condition is less stringent than Eq.(22c) (unity term) at low Ra and more stringent than Eq.(22a) at high Ra. We choose to drop the condition Eq.(22a) and to merge Eq.(22c) (unity term only) and Eq.(26a) into :

∆ h T 0 × P rRe 4N u 2 + P r 2 Re 3 Raθ ⋆ ∂t ≪ 1 (27)
Up to a numerical factor (1 -γ -1 ) of order one, the second thermal condition Eq.( 26b) is equivalent to the criterion (24). We now consider one by one each term in the brackets of Eq.(26c). The unity term gives a condition which is less stringent than Eq.( 27). The αT 0 term gives a condition less stringent than Eq.(26b) and can be ignored. The ∆ h /∆θ ⋆ rms term gives a condition comparable, or more stringent for high Ra, than Eq.( 23). We therefore keep the former one and re-write it alone as

∆ 2 h T 0 ∆ × P rRe 2θ ⋆ ∂t ≪ 1 (28)

D. Summary of the criteria and discussion

The set of the most stringent criteria for the applicability of Boussinesq approximation leaves us with only four sufficient conditions Eq.(20a), ( 24), ( 27) and (28), which are recalled below, up to a factor 2.

α∆ ≪ 1 (29a)

α∆ h × P rRe 2 Raθ ⋆ rms ≪ 1 (29b) ∆ h T 0 × P r 2 Re 3 Raθ ⋆ ∂t + P rRe 4N u 2 ≪ 1 (29c) ∆ 2 h T 0 ∆ × P rRe 2θ ⋆ ∂t ≪ 1 (29d)
and we recall that :

α∆ h = α 2 gT 0 h/c p = ρ 0 gχh(1 -γ -1 ) ∼ ρ 0 gχh
These 4 criteria are a first outcome of this paper. Although the choices of these specific criteria is not free from some arbitrariness, the order of magnitude of the LHS of these inequalities should be fairly robust to these choices.

Physically, the first criterion is required for incompressibility in the boundary layer. The other criteria are more difficult to relate to a specific physical effect. For example, the complete heat equation has a term accounting for the thermo-mechanical energy coupling which directly results from stratification [START_REF]This thermo-mechanical effect is already significantly compensated in the Rayleigh-Bénard geometry by the use of potential temperature[END_REF]. This term -which appears on the LHS of Eq.(12c)-gives contributions proportional to α∆ h , ∆ h /T 0 and ∆ 2 h /T 0 ∆. The fourth criterion above accounts for the ∆ 2 h /T 0 ∆ contribution alone. The third criterion above includes the ∆ h /T 0 contribution but also another condition necessary to neglect the viscous dissipation with respect to advected heat flux. The second criterion above allows to neglect variations of density due to the pressure variations, versus those due to temperature variations (Eq. (24). It also allows to neglect the cooling/heating associated with the pressure variations experienced by a fluid particle in the heat transport equation (Eq.(22b) and Eq.(26b)). These 4 criteria don't allow to assess a-priori the fulfilment of the Boussinesq approximation for given Ra and P r since it requires some knowledge of Re, θ ⋆ ∂t and θ ⋆ rms . A further step in the analysis consists in replacing these quantities by fits versus Ra and P r. For example, using the fits presented previously, the 4 previous criteria becomes :

α∆ ≪ 1 (30a) 0.1Ra 0.13 (α∆ h ) ≪ 1 (30b)
Ra 0.1 P r -0.25 x 200 + 10Ra -0.17 P r 0.25

∆ h T 0 ≪ 1 (30c) 0.1Ra 0.1 P r 0.25 x ∆ 2 h T 0 ∆ ≪ 1 (30d)
where we rounded prefactors and exponents to recall that the above inequalities are derived from order-ofmagnitude fits. We recall that the order of magnitude of x is taken as unity for P r of order 1-10. These equations highlight parameters having a direct influence on the validity of the Boussinesq approximation. For example, the dependence versus Ra is very weak. The equations Eq. (30) indicates that high Ra are weakly unfavourable but the uncertainty in the fits' exponents prevents from drawing precise conclusions on this weak dependence. For given Ra and P r, we found general features. Tall cells are unfavourable because of stratification (linear and quadratic dependence versus h coming from ∆ h ). Close proximity to the critical point is unfavourable too since α has a strong divergence in a region where ∆ h stays approximately constant. Finally, large temperature differences ∆ have opposite influence on the criterion Eq.(30a) (unfavourable) and on the criterion Eq.(30d) (favourable). In practice all these parameters can't all be varied independently.

Numerical applications are now proposed for a few typical laboratory Rayleigh-Bénard cells. Let's first consider a water cell of height h = 1 m, mean temperature T 0 = 50 • C and driven with ∆ = 30 • C (near atmospheric pressure). Such cells are being studied -for example-in physics labs in Lyon, Hong-Kong and Santa Barbara (for example see [54,55,56]). We found that the criteria Eq.( 30) are all fulfilled. The most stringent one is by far Eq.(30a), which is nevertheless still valid with 2 decades separation between the two sides of the inequality. The three other criteria are fulfilled with at least 5 decades of separation.

In h = 8 m air cell at atmospheric pressure, such as the one operated in Ilmenau (for example see [57]), we also find that the first criterion Eq.(30a) is the most stringent. For ∆ = 50 • C and T 0 = 50 • C, we find a ration 6 separation between the two sides of this inequality. The other criteria are fulfilled with 2.5 or more decades of separation.

Although atmospheric convection differs from RB convection of dry air, it is interesting the examine the criteria Eq.( 30) for h = 10 km, ∆ = 50 C and T 0 = 0 • C. Using dry air properties, we find that Eq.(30b) is violated by more than one decade and than the LHS of Eq.(30c) is close to 1, indicating that the Boussinesq approximation is not valid in this academic case. Surely these numbers should be taken with much caution because they rely on the extrapolation up to Ra ≃ 2.10 21 of fits derived at much lower Ra and probably in a different turbulent regime.

The case of Helium cells, such as those being operated in Grenoble and Trieste cannot be addressed with the same general procedure as above because, from one run to another, the temperature and the fluid density are varied and the parameter α does change accordingly (while ∆ h /h remains close to 2.10 -3 -3.10 -3 Km -1 ). Besides, the modelling presented in this paper ignores the possibility of occurrence of the Kraichnan's ultimate regime of convection [32], while the interpretation of some high Ra data is done according to this regime. These words of caution written, we can still get guidelines from the derived criteria. We numerically find that the two first criteria Eq.(30a) and Eq.(30b) are the most stringent ones, while the two others are fulfilled by at least 3 orders of magnitude for cells of height 20 cm and 1 m Therefore, we focus below on the validity of the second criterion for experimental data fulfilling the first one.

Data from [24,25,26] have been corrected from the side wall effect [58,59,60] and filtered to reject those with α∆ > 0.2, that is a 10% maximum relative variation of the density in each boundary layer (first criterion). All cells are cylindrical and their height and aspect ratio Γ = diameter/height are reported in the legend of figure 1.

The top plot shows the LHS of the criterion Eq.(30b) versus Ra. We recall here that the fit of θ ⋆ rms used to derive Eq.(30b) from Eq.(29b) is smooth over nearly 6 decades, up to nearly Ra = 10 15 and that we assume that its extrapolation over one extra decade (up to Ra = 10 16 ) provides a correct order of magnitude. From the plot, we find that this non-Boussinesq parameter reaches values larger than 0.2 for the highest Ra (Ra ∼ 10 16 ), which may no longuer be considered as significantly smaller than unity. This suggests that Boussinesq approximation could be violated for such ultra high Ra.

Besides, this shows that increasing the He cell height h to reach higher Ra -in Boussinesq conditions-is at best efficient like h 2 and not h 3 . Indeed, if we neglect the weak Ra dependence of Eq.(30b) and increase by a factor -say-10 the cell height, α should be 10 times smaller for α∆ h to stay constant (fulfilment of Eq.(30b)). On the other hand, ∆ can be made 10 times larger (fulfilment of Eq.(29a)). In temperature-density area where these high Ra data are obtained, a 10 fold decrease of α is roughly correlated with a 10 fold increase of the thermal diffusivity but with little change in the kinematic viscosity. All parameters considered, the maximum "Boussinesq" Ra will only be 100 times larger. One the other hand, this implies that such tall cells could be interesting models for atmospheric convection. Indeed, according to the estimate presented above we expect a similar violation of the Boussinesq approximation in the atmosphere.

The bottom plot of this figure illustrates that the compensated Nusselt N uRa -1/3 has different trends for Ra > 10 12 as mentioned earlier in the text. Regarding this issue, the top plot also shows that the non-Boussinesq parameter is below 10 -2 for Ra ∼ 10 13 and -more important-that it is approximately the same for datasets with different N u(Ra) trends (the vertical line across the plots illustrates this point). This shows that non-Boussinesq deviations of the sort encompassed by the present analysis cannot explain the discrepancy of heat transfer efficiency among He cells for 10 12 < Ra < 10 14 [23,24,25,26,27,28,31].

E. Conclusion and Perspectives

We derived a set of 4 criteria (Eq.( 29)) for assessing the applicability of Boussinesq approximation in a fluid with constant properties at high Rayleigh number Ra and intermediate P r. The fourth criterion, scaling quadratically with stratification, was found to be always less stringent than the 3 others in our numerical estimations.

The first criterion α∆ ≪ 1 simply reflects the incompressibility condition in the boundary layers. The 3 others are indirectly associated with stratification effects. They are written as functions of the adiabatic gradient but the thermodynamics relation Eq.( 6) allows to translate them into equivalent density stratification criteria. Contrary to a widespread idea, tall cells are not necessarily more Boussinesq than small ones for given Ra and P r, because of these stratification effects.

Using fits derived from the RB literature, we found that the first criterion α∆ ≪ 1 is the most relevant one for assessing the "Boussinesqness" of typical laboratory experiments conducted with water and air. In geophysical flows for which stratification occurs on very large scales, the second and third criteria could become more relevant than the first one.

In cryogenic helium, in addition to the first criterion, the validity of second one (Eq.( 29b)) deserves to be carefully examined at ultra high Ra, above 10 16 typically. This point and its consequences on very high Ra experiments is discussed in the text. At such high Ra, this second criterion is a sufficient condition for density fluctuations in the bulk to be more correlated to temperature fluctuations than to pressure fluctuations (Eq.( 24)) and to neglect, in the heat equation, a thermal heat source (or sink) term due to pressure-fluctuation-induced compression (or expansion) (Eq.(26b)).

At lower Rayleigh numbers (typically 10 12 < Ra < 10 14 ), we find that the apparent incompatibility of heat transfer measurements reported in literature cannot be explained by such non-Boussinesq effects.

In the Boussinesq approximation, the equations of the fluid are only parametrised by Ra and P r. The present work shows that a weaker approximation of the equations -preserving some terms proportional to the parameters α∆ and/or α∆ h (or equivalently χρgh)-would enable to account for first deviations from the Boussinesq limit of turbulent Rayleigh-Bénard cells filled with a fluid with constant properties. As a perspective, a DNS simulation conducted with such generalised equations at moderate turbulent Rayleigh numbers (say 10 8 < Ra < 10 10 ) but with large enough values of the non-Boussinesq parameters α∆ or α∆ h would provide valuable quantitative information on deviations from the Boussinesq limit. In particular, simulation with α∆ h of order unity can still be performed assuming incompressibility. Such simulations would give a precise quantitative meaning on the ≪ symbol used in the proposed Boussinesq criteria.

FIG. 1 :

 1 FIG.1: Data from three cryogenic He experiments[24,25,26], restricted to those fulfilling α∆ < 0.2. The top plot shows the parameter appearing on the left of the Boussinesq criterion Eq.(30b). The Boussinesq condition is better fulfilled for the values of the parameter smaller than one (the Chavanne et al. data near Ra ∼ 10 12 sitting below 10 -4 have been obtained in liquid helium). The bottom plot shows the compensated Nusselt N uRa -1/3 versus Ra. We note that, for Ra ∼ 10 13 (vertical line), the datasets already sits on different N u(Ra) trends while some of their non Boussinesq parameters, displayed on the top plot, roughly overlap.

Acknowledgments

We thank R. Verzicco for sharing DNS data, which were used to validate a scaling law for pressure fluctua-tions, Q. Zhou for providing some physical properties of water, J. Niemela for sending data from his aspect ratio 0.5 experiment and my colleagues B. Chabaud and H. Hébral for proof-reading. We acknowledge discussion with M. Pons, P. Carles and more especially with B. Castaing. This work was initiated thanks to the Région Rhône-Alpes contract 301491302.


Babinet-bilayered geometric phase optical elements

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

Babinet-bilayered geometric phase optical elements

ETIENNE BRASSELET

Univ. Bordeaux, CNRS, LOMA, UMR 5798, F-33400 Talence, France (etienne.brasselet@u-bordeaux.fr) Geometric phase optical elements based on structured anisotropy are widely used for phase shaping via their orientational degree of freedom. To date, amplitude shaping via space-variant retardance is much less investigated, a practical reason being that the spin-orbit interaction of light couples retardance with the dynamic part of the optical phase. Inspired by the complementary diffractive elements associated with Babinet's principle, a bilayered subwavelength grating design is proposed in order to cancel out the spatial modulation of the dynamic phase usually associated with space-variant birefringent phase retardation. This concept is illustrated in the framework of single-mode Laguerre Gauss -beam shaping.

OCIS codes: (050.6624) Subwavelength structures; (140.3300) Laser beam shaping.

Structuring matter at the subwavelength scale is known to allow the design and manufacture of artificial materials with wave properties that are not possessed by the constituting bulk media taken apart. These materials are usually called metamaterials, a name derived from the Greek "meta" (beyond) emphasizing the access to novel properties and phenomena that seems limited so far only by our imagination. In electromagnetism, metamate rials evolved from radio frequencies to optics with substantial conceptual and engineering advances during the last two dec ades [1]. This concept went beyond electromagnetism, and has been addressed when thermodynamics, acoustics, elasticity, or quantum mechanics are at play [2]. Quite naturally, the process of progress has led to the development of dynamic versions of static metamaterials toward the creation of enhanced or novel functionalities [3,4] and now embarks the richness and subtleties of nonlinear physics as well [5,6].

Of course, advanced designs continue to be developed even in the linear case. For instance, in the context of spin controlled asymmetric orbital angular momentum states, it has been pro posed recently to encode independent dynamic and geometric phase profiles into single metallic [7] or dielectric [8] metasur faces. In practice, this is made by controlling locally the orien tation and aspect ratio of equally spaced subwavelength anisotropic nanostructures. Besides phase control, amplitude shaping is also accessible to metasurfaces by exploiting their vectorial nature. This requires the combination of the local con trol of the birefringent phase retardation with post polarization filtering, which has been exploited recently in the context of modal beam shaping [9]. In short, the optical characteristics of a passive metasurface are determined by the spatial distribu tion of three quantities: (i) the dynamic phase, Φ dyn , which is associated with the average optical path length; (ii) the geomet ric phase, Φ geom , which uniquely depends on the structure of the medium; and (iii) the birefringent phase retardation, Δ, which arises from anisotropic phase delay between the ordinary and extraordinary waves. From a general point of view, Φ dyn , Φ geom , and Δ can be independently adjusted by design. However, this requires a conception step that relies on optimiz ing these parameters via numerical simulations, which prevents from an intuitive approach to conceive novel photonic functionalities.

Here, inspired by the Babinet's principle, a bilayered design made of two complementary subwavelength gratings (see Fig. 1) is proposed with the aim at getting rid of the dynamic phase spatial modulation usually associated with amplitude modulation driven by space variant birefringence. This is done by multiplexing two independent structural degrees of freedom that independently act on amplitude and phase. Namely, we address the design of metasurfaces having a uniform Φ dyn while (i) the grating period and its filling factor define Δ, and (ii) the orientation angle, ψ, of the grating wavevector defines Φ geom 2σψ, where σ 1 is the helicity of the incident light. In order to gain insight, it is worth recalling that the behavior of a transparent optical metasurface in the paraxial regime can be grasped simply and analytically. Indeed, by con sidering an incident circularly polarized field described by the unit vector c σ x iσy∕ 2 p , where x, y is the plane of the metasurface, the output light field can be expressed as [10] E out ∝ e iΦ dyn cosΔ∕2c σ i sinΔ∕2e iΦ geom c -σ : (1)

The sought after Babinet metasurface is ideally character ized by an effective σ dependent contra circular complex transmittance, τ eff ∝ sinΔ∕2 exp2iσψ:

(2)

Qualitatively, the proposed Babinet bilayered design shown in Fig. 1 in the one dimensional case of a square waveform with subwavelength spatial period Λ can be understood from the fol lowing reasoning. Obviously, a slab of an isotropic dielectric is associated with a uniform dynamic phase. In contrast, once peri odically structured, it behaves as a uniaxial dielectric slab [11] with a dynamic phase that depends on the amount of the re moved material that is characterized by the filling factor F. The idea comes from the following intuition: the sequence of two complementary subwavelength gratings obtained by split ting an isotropic dielectric slab is thought as an effective subwave length structured slab without matter removal overall. As such, one could expect a constant dynamic phase independently of F , while the form birefringence remains strongly dependent on F .

The above reasoning is investigated quantitatively in the framework of the second order effective theory [12], by con sidering diamond material in air as a case study, whose refractive indices are, respectively, taken as N 2.425 and 1 at λ 532 nm wavelength. More precisely, a subwavelength grating illuminated at normal incidence and characterized with a filling factor F and a period Λ behaves as a uniaxial medium with refractive indices along a direction parallel (∥) and perpendicular (⊥) to the grating wavevector that are given by

n ∥ n ∥,0 1 π 2 Λ 2 3λ 2 F 2 F 2 n 4 ∥,0 n 2 ⊥,0 1 1 N 2 2 1∕2 , (3)
n ⊥ n ⊥,0 1 π 2 Λ 2 3λ 2 F 2 F 2 N 2 1 2 1∕2 , ( 4
)
where λ is the wavelength, F 1 F , n ⊥,0 N 2 F F 1∕2 , and n ∥,0 N F N 2 F -1∕2 . Noting that the dynamic phase of the bilayer is Φ dyn 2kL ñ, where k 2π∕λ and L is the thickness of each layer and ñ n n∕2 with n and n being the average re fractive index of the first subwavelength grating and its comple ment, the dynamic phase behavior as a function of the structural characteristics is assessed via the study of ñ. This is reported in Fig. 2(a) where ñ is plotted in the plane of param eters Λ, F . Here, the chosen range 0 < N Λ∕λ < 1 refers to the absence of diffraction orders other than the zeroth order one, which practically defines a subwavelength grating. Clearly, there is an optimal value for the grating period value that ensures a minimal dependence of ñ on F , as illustrated in Fig. 2(b). Namely, as shown in the inset of that figure, N Λ∕λ 0.789 minimizes the standard deviation of ñF , which is the value used further in this work. The optical proper ties of the individual subwavelength gratings constituting the optimized bilayered element are summarized in Fig. 3, where the effective extraordinary (n ∥ and n∥ ), ordinary (n ⊥ and n⊥ ), and average refractive (n and n) indices are plotted as a function of F . These two plots emphasize the compensation mechanism leading to an effective average refractive index of the bilayered optical element that varies much less with F than the average refractive index of any of the two gratings taken separately. On the other hand, the birefringence of the complementary gratings have similar behavior, hence adding up.

Said differently, the filling factor can be used to modulate the birefringent phase retardation of the bilayer, Δ 4πd ñL∕λ with d ñ ñ∥ ñ⊥ and ñ∥,⊥ n ∥,⊥ n∥,⊥ ∕2, without much variation of the dynamic phase. This is summarized in Fig. 4 where the effective average refractive index and bire fringence are plotted as a function of F. As shown in Fig. 4(a), ñ exhibits slight modulations near F → 0 and F → 1, which implies that the optical element suffers from residual modula tion of the dynamic phase if one exploits the full range for F. Nevertheless, considering a situation that corresponds to a birefringent phase retardation Δ π at F 0.5, there is a substantial decrease of the magnitude of the dynamic phase modulation in the case of the Babinet bilayer compared to the case of a single layer. This is illustrated in Fig. 5, where the variation of the dynamic phase in the range 0 < F < 0.5 is shown in both cases, imposing a fixed maximal birefringent phase retardation at F 0.5. The dynamic phase profile is flat tened by a factor of ≃42 when considering the full range 0 < F < 0.5, and this factor increases up to ≃395 by restrict ing to the range 0.2 < F < 0.5, hence limiting the drawbacks associated with the resolution capabilities of nanofabrication tools as F → 0. However, in order to obtain a full modulation of τ eff between 0 and 1, one has to increase L in order to fulfill the condition Δ 2π (i.e., τ eff 0) at F 0.5, which gives a thickness L ≃ 450 nm and Δ π (i.e., τ eff 1) at F ≃ 0.17 with the present parameters. Since the condition N Λ∕λ 0.789 corresponds to Λ ≃ 170 nm, the latter conditions in volve the fabrication of subwavelength gratings with aspect ratio L∕F Λ up to 15, which is accessible to nanofabrication techniques [13].

Next, the interest of the proposed concept is illustrated in the framework of beam shaping by considering the conception of geometric phase Laguerre Gaussian (LG) modal beam shapers.

LG beams refer to a complete orthogonal basis for the scalar para xial Helmholtz equation [START_REF] Siegman | Lasers[END_REF], each element of the basis being univocally associated to an azimuthal index, l , and a radial index, p. Many strategies have been developed toward the generation of LG like beams using either extra cavity or intra cavity schemes, as reviewed in [9]. The use of geometric phase optical elements is especially interesting in that the helicity dependent geometric phase offers a route to switch between a mode l , p and a mode l , p by flipping the incident helicity of light from σ to σ, provided that ψ σl ∕2ϕ.

Since the first LG like experimental demonstrations in the case p 0 in the early 2000s [START_REF] Biener | [END_REF] and p ≥ 1 only a few years ago [16,17], the creation of geometric phase LG modal con verters as such has not yet been reported. Recently, a single layer theoretical proposal combining dynamic and geometric phase modulations at fixed n ⊥ has been reported [9]. However, the corresponding design breaks the right/left sym metry and thus works for only a single value of the helicity. A quasi modal version recovering the latter broken symmetry has been discussed on experimental grounds [18], though at the expense of substantial contra circular transmission losses.

Here, we show that the use of space variant Babinet bilayered subwavelength gratings allows to create polarization controlled LG modal beam shapers. This is discussed in the representative case p 0; hence, the index p is omitted in what follows, and we note that the generalization to high order radial modes poses no particular problem except more cumbersome calculations. At first, we recall that a LG l field associated with the waist radius w is expressed in its focal plane as

E w l ∝ r w jl j exp r 2 w 2 ilϕ : (5)
Then we deduce from Eq. ( 2), which describes an ideal Babinet bilayer optical element, that the use of a birefringence phase retardation radial profile Δ l r satisfying the condition

sinΔ l r∕2 r jl j exp r 2 w 2 r 2 w 2 in max r h r jl j exp r 2 w 2 r 2 w 2 in i , (6)
combined with an optical axis orientation angle

ψ l ϕ l ∕2ϕ, (7)
ensures the generation of a LG σl mode with waist w by placing the optical element in the focal plane (or, in practice, by using a collimated beam) of a normally incident circularly polarized Gaussian beam with waist w in , provided that w < w in . The efficiency of the device (i.e., the fraction of the incident power that is converted into the desired mode power) is 4∕w 2 in R ∞ 0 sin 2 Δ l r∕2 exp 2r 2 ∕w 2 in rdr and depends on w∕w in . Its optimal value is reached for w opt w in 1 jl j -1∕2 . In this case, introducing ρ r∕w in , Eq. ( 6) simplifies to sinΔ opt l ρ∕2 2 jl j∕2 ρ jl j exp jl jρ 2 jl j∕2:

The modal figure of merit of a space variant Babinet bilayer design given by Eqs. ( 7) and ( 8) is assessed quantitatively by evaluating the optimal power fraction C σl of the σ polarized field component at the output of the element that corresponds to a LG σl mode. Namely, where E out is given by Eq. ( 1), and the asterisk denotes complex conjugation, paying attention to the fact that LG modal decom position is not unique [19,20]. By considering the retardance range 0 < Δ < π, we obtain a high degree of modality. Namely, C l > 0.998 for l 1,2, 3, which is compared to two other situations. First, the case of single subwavelength gra tings characterized by the same structural features as the Babinet bilayers regarding both the optical axis orientation angle [Eq. ( 7)] and the birefringence phase retardation [Eq. ( 8)], to which we refer as the "non twisted single layer" case. Second, the case discussed in Ref. [9], to which we refer as the "twisted single layer" case, according to ψ l r, ϕ l ϕ Φ dyn r∕2. In that case, the bire fringence phase retardation is also given by Eq. ( 8), while the op tical axis orientation angle is designed to cancel out the dynamic phase for one of the two helicity states (namely, C l 1) at the expense of doubling it for the other (namely, C -l < 1).

C σl w R ∞ 0 E w l E out • c -σ rdr 2 R ∞ 0 jE w l j 2 rdr R ∞ 0 jE out • c -σ j 2 rdr , (9)
The results are summarized in Fig. 6 and Table 1, both emphasiz ing the superior modal performances of Babinet bilayers.

In the field of optical information processing, the spin re versal symmetry of the proposed approach allows considering the extension of polarization controllled high order Poincaré spheres to the case of LG modes with given azimuthal and radial indices, l and p. This solves a serious limitation of present approaches restricted so far to an infinite superposition of radial modes at fixed l , based on non modal geometric phase elements associated with ψ l ∕2ϕ and Δ π [21].

In practice, there is no need for stringent subwavelength precision relative orientational alignment of the two subwave length gratings, as we are dealing with effective media. In addition, a monolithic design consisting of placing the two structured layers at both sides of a substrate can be considered, whose feasibility is supported by the previously reported fab rication of double sided geometric phase lenses [22]. Still, one should keep in mind that the grating depth and the nature of the surrounding media play important roles in the effective properties of the structured material [23]. In addition, Fresnel reflections may also be taken into account, at least for two main reasons: (i) the modification of the incident polarization as light enters into the structures due to optical anisotropy and (ii) cav ity effects associated with the presence of two structures. These possible drawbacks are evaluated at F 0.5. This gives for the case (i), the spin flipping of ≲10% of the incident photons for incident circular polarization, and for the case (ii), a typical Fabry Pérot finesse ≲1. We conclude that although an attempt to implement such a design is not excluded in principle, quantitative attention must be kept in mind.

Finally, we note that it is not the first time that Babinet's prin ciple has inspired the physics of metasurfaces and metamaterials, a pioneering example being the introduction of complementary split ring resonators acting as electric point dipoles with negative polarizability [24], which was followed by a huge number of de velopments. Present work therefore illustrates how basic optical concepts can bring attractive photonic opportunities to a field that has much progressed during the last two decades [25], in particular when the orbital angular momentum is at play. a This case refers to Ref. [9], and the two values refer to σ 1.

b These cases are independent on σ.

Fig. 1 .

 1 Fig. 1. (a) Slab of an isotropic dielectric with thickness L. (b) Proposed Babinet bilayered design that consists of two comple mentary subwavelength gratings made by splitting the slab shown in panel (a). By construction, the two complementary gratings have identical periods, Λ, and complementary filling factors, F and F 1 F .

Fig. 2 .

 2 Fig. 2. (a) Effective average refractive index ñ of the Babinet bilay ered subwavelength grating versus the reduced spatial period N Λ∕λ and the filling factor F. (b) Maximal ( ñmax max ñF ), minimal ( ñmin min ñF ), and average (h ñi F ) effective refractive indices ver sus N Λ∕λ. Inset: standard deviation of ñF as a function of the reduced period whose value N Λ∕λ 0.789 (see arrow) defines the condition ensuring a minimal dependence of ñ on F .

Fig. 3 .

 3 Fig. 3. (a) Extraordinary (n ∥ ), ordinary (n ⊥ ), and average (n) refrac tive indices of a subwavelength grating as a function of the filling factor F . (b) Same as in (a) for the complementary subwavelength grating. The simulations are made for N Λ∕λ 0.789, which optimizes the dynamic phase flatness.

Fig. 4 .

 4 Fig. 4. (a) Extraordinary ( ñ∥ ), ordinary ( ñ⊥ ), and average ( ñ) effective refractive indices of the optimized Babinet bilayered subwavelength grating that corresponds to N Λ∕λ 0.789 according to the inset of Fig. 2(b). (b) Corresponding effective birefringence, d ñ ñ∥ ñ⊥ .

Fig. 5 .

 5 Fig. 5. Normalized variation of the dynamic phase as a function of the filling factor for a single subwavelength grating (black curve) and its comparison with the optimized Babinet bilayer (red curve) for a given maximal retardance at F 0.5.

Fig. 6 .

 6 Fig. 6. Modal coefficients C l versus the waist w of the sought afterLG σl mode for three modal designs and l 1,2, 3.

Table 1 .

 1 Optimal Modal Coefficients max w C σl w for Three Different Designs and l 1,2, 3, which Corresponds to the Maximal Values Associated with the Plots Shown in Fig.6

	Grating Type	' 1	' 2	' 3
	Twisted monolayer a	(1,0.082)	(1,0.274)	(1,0.520)
	Non twisted monolayer b	0.454	0.555	0.744
	Babinet bilayer b	0.999	0.998	0.998