Mondher Bouteraa

Chérif Nouar
email: cherif.nouar@univ-lorraine.fr

Weakly nonloinear analysis of Rayleigh-Bénard convection in a non-Newtonian fluid between plates of finite conductivity: Influence of shear-thinning effects

Keywords: bifurcation, Rayleigh-Bénard convection, shear-thinning fluid

Finite-amplitude thermal convection in a shear-thinning fluid layer between two horizontal plates of finite thermal conductivity is considered. Weakly nonlinear analysis is adopted as a first approach to investigate non linear effects. The rheological behavior of the fluid is described by the Carreau model. As a first step, the critical conditions for the onset of convection are computed as a function of the ratio ξ of the thermal conductivity of the plates to the thermal conductivity of the fluid.

. For a significant deviation from the critical conditions, nonlinear convection terms and nonlinear viscous terms become stronger leading to a further diminution of ξ c . The dependency of the heat transfer on ξ and the rheological parameters is reported. It is consistent with the maximum heat transfer principle. Finally, the flow structure and the viscosity field are represented for weakly and highly conducting plates.

I. INTRODUCTION

Recently, a weakly nonlinear stability analysis of thermal convection in a layer of a non-Newtonian fluid between two horizontal plates heated from below was considered by Bouteraa et al. [START_REF] Bouteraa | Weakly nonlinear analysis of Rayleigh-Bénard convection in shear-thinning fluids: nature of the bifurcation and pattern selection[END_REF]. Assuming the fluid purely viscous and shear-thinning, Bouteraa et al. [START_REF] Bouteraa | Weakly nonlinear analysis of Rayleigh-Bénard convection in shear-thinning fluids: nature of the bifurcation and pattern selection[END_REF] studied the influence of shear-thinning effects on the nature of the bifurcation, the convection intensity and the pattern selection. The possibility of wall slip was taken into account by using Navier's slip law at top and bottom walls. It was shown that the bifurcation is supercritical for moderately shear-thinning effects and becomes subcritical for strongly shear-thinning effects. The critical value of the degree of shear-thinning α c defined by Eq. ( 14), above which the bifurcation becomes subcritical is determined as a function of a dimensionless slip length parameter. Bouteraa et al. [START_REF] Bouteraa | Weakly nonlinear analysis of Rayleigh-Bénard convection in shear-thinning fluids: nature of the bifurcation and pattern selection[END_REF] demonstrated that near the threshold of the convection, only rolls are stable and this stability is reinforced by the shear-thinning behavior. It is shown that shear-thinning effects increase the interaction between sets of rolls that constitute square and hexagon patterns. In the above analysis, it was assumed that the plates are of infinite thermal conductivity kp and heat capacity.

In experimental situations, however, the thermal conductivity of the plates is finite. Furthermore, in some situations, the plates are much poorer conductors than the fluid. For instance, in LeGal & Croquette [START_REF] Legal | Appearance of a square pattern in a Rayleigh-Bénard experiment[END_REF], the plates are made of Plexiglas and the fluid is water.

The ratio ξ of the thermal conductivity of the plates kp to the thermal conductivity of the fluid k is ξ = 0.4. In Kebiche [START_REF] Kebiche | Etude expérimentale de l'instabilité de Rayleigh-Bénard dans les fluides non-Newtoniens[END_REF] and Kebiche et al. [START_REF] Kebiche | Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluids[END_REF], the horizontal plates are made of Polycarbonate and the fluid is an aqueous solution of CarboxyMethylCellulose or Carbopol, the ratio ξ is estimated to ξ = 0.25. In these situations, as well as those encountered in some engineering convection problems and in geophysics, the assumption that the plates are held at fixed and uniform temperatures loses its validity. Indeed, temperature disturbance penetrates the wall. More general boundary conditions have to be used. The temperature and the heat flux obey a continuity equation at the boundaries.

For a Newtonian fluid, the effect of imperfect conducting plates on the Rayleigh-Bénard convection was investigated by several authors. First, in the linear theory frame, Sparrow et al. [START_REF] Sparrow | Thermal instability in a horizontal fluid layer: effect of boundary conditions and a non-linear temperature profile[END_REF] and Hurle et al. [START_REF] Hurle | On the solution of the Bénard problem with boundaries of finite conductivity[END_REF] (see also Riahi [START_REF] Riahi | Nonlinear convection in a porous layer with finite conducting boundaries[END_REF], Clever & Busse [START_REF] Clever | Convection in a layer heated from below with a nearly insulating boundary[END_REF] and Holmedal Rayleigh number Ra c and the critical wave number k c , vary continuously from 1708 to 720 and from 3.11 to 0, when the ratio of thermal conductivities decreases from infinity to zero.

According to Cerisier et al. [START_REF] Cerisier | Thermal influence of boundaries on the onset of Rayleigh-Bénard convection[END_REF] a temperature fluctuation occurring in the liquid close to a nearly insulating plate persists and distorts the temperature distribution. This distribution can led to instability of the fluid layer. As a consequence, the temperature gradient is small and the fluid organizes in a pattern with a small wavenumber.

Exploiting the fact that for nearly insulating walls (ξ << 1), the horizontal scale of convection is much larger than the depth of the fluid, Busse & Riahi [START_REF] Busse | Nonlinear convection in a layer with nearly insulating boundaries[END_REF] considered weakly nonlinear three-dimensional solutions in the case of infinitely thick plates. They showed that near the onset, square convection cells are the stable planform, in contrast with two-dimensional rolls which are the only stable convection pattern in a symmetrical situation with isothermal boundary conditions. This result was confirmed and extended to fully nonlinear convection by Proctor [START_REF] Proctor | Planform selection by finite amplitude thermal convection between poorly conducting slabs[END_REF] using a "shallow water theory" adapted for the Rayleigh-Bénard convection by Chapman & Proctor [START_REF] Chapman | Nonlinear Rayleigh-Bénard convection between poorly conducting boundaries[END_REF]. The thickness of the plates was assumed finite. Jenkins & Proctor [START_REF] Jenkins | The transition from roll to square-cell solutions in Rayleigh-Bénard convection[END_REF] considered three-dimensional finite-amplitude thermal convection with finite thickness and finite thermal conductivity of the bounding plates. They determined the critical value ξ c of the thermal conductivities ratio at which the preferred planform changes from square cell (ξ < ξ c ) to roll (ξ > ξ c ), as function of the Prandlt number. When the thickness of the plates is of the same order as that of the fluid layer and for P r ≥ 10, rolls should be observed when ξ > 1 and squares when ξ < 1. This is in agreement with LeGal & Croquette's experiments [START_REF] Legal | Appearance of a square pattern in a Rayleigh-Bénard experiment[END_REF].

The objective of the present paper is twofold. Firstly, to determine the critical value of the shear-thinning degree above which the bifurcation becomes subcritical as a function of ξ. Secondly, to study the influence of shear-thinning effects on ξ c . We hope that our findings will shed new light on the interpretation of the results obtained by Kebiche [START_REF] Kebiche | Etude expérimentale de l'instabilité de Rayleigh-Bénard dans les fluides non-Newtoniens[END_REF] and Kebiche et al. [START_REF] Kebiche | Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluids[END_REF], although the fluid used in [START_REF] Kebiche | Experimental investigation of the Rayleigh-Bénard convection in a yield stress fluids[END_REF] is not only shear-thinning but has also a yield stress.

The paper is organized as follows. In section II, the mathematical formulation of the problem is presented. In section III, the linear stability analysis for the onset of convective flow is reinvestigated. The critical Rayleigh number (Ra c ) and wavenumber (k c ) are nonlinear stability analysis. The results are discussed in section V. The critical value of the shear-thinning degree above which the bifurcation becomes subcritical is determined as a function of ξ. The pattern selection near the threshold of convection is investigated in terms of α and ξ. For a significant deviation from the critical conditions, higher-order solutions are computed in section VI. Section VII provides information on the flow structure and the heat transfer. A concluding discussion is given in the last section of the paper.

II. PROBLEM FORMULATION

A. General equations and parameters

We consider a layer of shear-thinning fluid of depth d confined between two horizontal plates that are infinite in extent and which have a thickness Λ d, where Λ is of order unity. The outer surface of the bottom and top plates are kept at constant temperatures respectively T0 + ∆ T /2 and T0 -∆ T /2, with ∆ T > 0. The fluid has density ρ, thermal conductivity k, thermal coefficient expansion (at constant pressure) β and viscosity μ0 at zero shear rate. The thermal conductivity of the slabs is kp . Because of the thermal expansion, the temperature difference between the two plates, induces a vertical density stratification.

Heavy cold fluid is above a light warm fluid. For small ∆ T , the fluid remains at rest and the heat is transferred by conduction. In the fluid, 0 < ẑ < d, the hydrostatic solution and the temperature profile are:

d P dẑ = -ρ 0 ĝ 1 -β T -T0 and Tcond = T0 + ∆ T 1 + 2Λ/ξ 1 2 - ẑd , (1)
where, ĝ is the acceleration due to gravity. Here, the z-axis is directed upwards, with the origin located at the bottom plate. The reference temperature T0 is the temperature in the middle of the fluid layer and ρ0 is the fluid density at T0 . The temperature difference between the top and the bottom of the fluid layer is ∆ Tf = ∆ T /(1 + 2Λ/ξ).

The temperature profile in the top and bottom plates are:

Tcond = T0 + ∆ T ξ + 2Λ 1 - 1 2 ξ - ẑd , d ≤ ẑ ≤ (1 + Λ) d (2)
and

Tcond = T0 + ∆ T ξ + 2Λ 1 2 ξ - ẑd , -Λ d ≤ ẑ ≤ 0. ( 3
)
When ξ is low, a large part of ∆ T occurs across the plates, and remains only a small part ∆ Tf of ∆ T , acting as the driving force for the convection. When ∆T f , exceeds a critical value, the buoyancy force overcomes, the dissipative effects, i.e., viscous and heat diffusion, the convection sets in and a so-called convective patterns emerge. The stability of the hydrostatic solution is considered by introducing temperature and pressure perturbation as well as a fluid motion. Boussinesq approximation is adopted, i.e., the temperature dependence of the fluid properties can be neglected except for the temperature induced density difference in the buoyant force that drives the flow. The heat production due to viscosity is neglected.

Here and in what follows, the quantities with hat (.) are dimensional. Distances are scaled with d, velocity with κ/ d, where κ is the thermal diffusivity of the fluid, time with d2 /κ (characteristic time scale of thermal diffusion), temperature with ∆ Tf , pressure and stresses with κμ 0 / d2 . Using these scales, the dimensionless perturbation equations read:

∇ • u = 0, (4) 1
P r ∂u ∂t + u • ∇u = -∇p + Ra θ e z + ∇ • τ , (5)
∂θ ∂t + u • ∇θ = u • e z + ∇ 2 θ , (6)
in the fluid, and

∂ θ ∂t = κp κ ∇ 2 θ , (7)
in the bounding slabs. Here, e z denotes the unit vector in the vertical direction, u(x, t) = ue x + ve y + we z is the fluid velocity and p(x, t) and θ(x, t) represent the pressure and temperature deviations from their values in the conduction state. The temperature perturbation in the slabs is denoted θ(x, t). The position vector x has components x, y, z. The Rayleigh number Ra and the Prandtl number P r are defined by :

Ra = ρ0 ĝ β∆ Tf d3 κ μ0 ; P r = μ0 ρ0 κ . (8)
Generally, for non-Newtonian fluids, P r >> 1, i.e., the viscous diffusion time is shorter than the thermal diffusion time.

The fluid is assumed to be purely viscous and shear-thinning. The viscous stress-tensor

τ = µ (Γ) γ with γ = ∇u + (∇u) T (9)
the rate-of-strain tensor, of second invariant

Γ = 1 2 γij γij . ( 10
)
The Carreau model is given by

μ -μ∞ μ0 -μ∞ = 1 + λ2 Γ nc-1 2 , ( 11
)
with μ0 and μ∞ the viscosities at low and high shear rate, (n c < 1) the shear-thinning index, λ the characteristic time of the fluid. The location of the transition from the Newtonian plateau to the shear-thinning regime is determined by λ, since 1/ λ defines the characteristic shear rate for the onset of shear-thinning. Increasing λ reduces the Newtonian plateau to lower shear rates. The infinite shear viscosity, μ∞ , is generally associated with a breakdown of the fluid, and is frequently significantly smaller (10 -3 to 10 -4 times smaller) than μ0 , see Bird et al. [START_REF] Bird | Dynamics of polymeric liquids[END_REF] and Tanner [START_REF] Tanner | Engineering rheology[END_REF]. The ratio μ∞ /μ 0 will be thus neglected in the following.

The dimensionless effective viscosity is then

µ = μ μ0 = 1 + λ 2 Γ nc-1 2 with λ = λ d2 /κ . ( 12
)
The Newtonian behavior, μ = μ0 , is obtained by setting n c = 1 or λ = 0.

For a small amplitude disturbance, the viscosity can be expanded about the hydrostatic solution,

µ = 1 + n c -1 2 λ 2 Γ + 1 2 n c -1 2
n c -3 2 λ 4 Γ 2 + ... (13)
At lowest nonlinear order, a relevant rheological parameter is the 'degree of shearthinning'

α = dµ dΓ Γ=0 = 1 -n c 2 λ 2 . ( 14
)
For the perturbation velocity u, the conditions imposed on the top and bottom plates are either no-slip (no-slip boundary conditions NSBC), which implies

u = 0 at z = 0, 1 (15)
or stress-free (stress-free boundary conditions SFBC), which implies

∂u ∂z = ∂v ∂z = w = 0 at z = 0, 1 (16)
For thermal boundary conditions, a constant and uniform temperature is assumed on the outer surface of each plate

θ = 0 at z = -Λ, 1 + Λ ( 17
)
The continuity condition for temperature and heat flux are

θ = θ at z = 0, 1, (18)
Dθ = ξD θ at z = 0, 1. (19)
Where D ≡ ∂ ∂z and ξ = kp k .

Remark:

The governing equations ( 5) and ( 6) with the constitutive equation ( 9) and the boundary conditions ( 15)-( 19) are reflection-symmetric about the midplane z = 1/2. The action of this so-called Boussinesq symmetry is

[u, v, w, θ, p] (t, x, y, z) → [u, v, -w, θ, p] (t, x, y, 1 -z).

D. Reduction: elimination of the pressure

In the momentum equations, the pressure field can be eliminated using the curl of Eq.

(5). We then take the curl of Eq. ( 5) one more time. Using the continuity equation, and

projecting onto e z , we get the following evolution equations for the vertical vorticity ζ and the vertical velocity w:

∂ζ ∂t + e z • ∇ × [(u • ∇) u] = P r ∆ζ + P r e z • ∇ × [∇ • (µ -1) γ] , (20)
∂∇ 2 w ∂t -e z • [∇ × ∇ × [(u.∇) u]] = P r ∆ 2 w + Ra P r ∇ 2 H θ - ( 21
) P r [∇ × ∇ × [∇ • (µ -1) γ]] • e z , ∂θ ∂t + (u • ∇) θ = w + ∇ 2 θ , (22)
∂ θ ∂t = κp κ ∇ 2 θ , (23)
where

ζ = ∂v ∂x - ∂u ∂y and ∇ 2 H = ∂ 2 ∂x 2 + ∂ 2 ∂y 2 .
From the continuity equation and the vertical vorticity definition, one can deduce the horizontal velocity components (u , v):

∇ 2 H u = - ∂ 2 w ∂x∂z - ∂ζ ∂y ; ∇ 2 H v = - ∂ 2 w ∂y∂z + ∂ζ ∂x . (24)
The boundary conditions for w are:

w = Dw = 0 at z = 0, 1 for NSBC, (25)
and

w = D 2 w = 0 at z = 0, 1 for SFBC. ( 26
)
For the temperature, the boundary conditions are :

θ = 0 at z = -Λ, 1 + Λ, (27)
θ = θ at z = 0, 1, (28)
Dθ = ξD θ at z = 0, 1. (29)
In the following, as in Chapman & Proctor [START_REF] Chapman | Nonlinear Rayleigh-Bénard convection between poorly conducting boundaries[END_REF], Proctor [START_REF] Proctor | Planform selection by finite amplitude thermal convection between poorly conducting slabs[END_REF] and Carriere et al. [START_REF] Carriere | Wavelength selection in Rayleigh-Bénard convection between horizontal boundaries of finite conductivity[END_REF], we will assume that κp /κ = ξ.

A. Critical conditions and critical modes

In the linear theory, u and θ are assumed infinitesimal. The nonlinear terms in ( 20)-( 22) can be neglected. We obtain :

1 P r ∂ζ ∂t = ∆ζ , (30) 1
P r ∂∆w ∂t = ∆ 2 w + Ra∆ H θ , (31)
∂θ ∂t = w + ∆θ , (32)
∂ θ ∂t = κp κ ∆ θ. ( 33
)
At this stage, no non-Newtonian effects enter the problem. The vertical vorticity decouples and obeys a diffusion equation and thus can be ignored in the linear theory. For equations ( 31)-( 33), we seek a normal mode solution

     w(x, y, z, t) θ(x, y, z, t) θ(x, y, z, t)      =      F 11 (z) G 11 (z) G11 (z)      f (x, y) exp (s t) , (34)
where s = s r +is i , is a complex number and f (x, y) satisfies the two-dimensional Helmoltz equation ∆ H f = -k 2 f . Here k is the norm of the horizontal wavenumber k. Substituting (34) into ( 31)-( 33), leads to the differential equations

s P r -1 D 2 -k 2 F 11 = -k 2 RaG 11 + D 2 -k 2 2 F 11 , (35)
s G 11 = F 11 + (D 2 -k 2 )G 11 , (36)
s G11 = κp κ (D 2 -k 2 ) G11 . ( 37
)
It may be shown easily that s is real. The boundary conditions are:

F 11 = DF 11 = 0 at z = 0, 1 For NSBC ( 38
)
F 11 = D 2 F 11 = 0 at z = 0, 1 For SFBC (39) G11 = 0 at z = -Λ, 1 + Λ, (40)
G 11 = G11 at z = 0, 1 (41)
DG 11 = ξD G11 at z = 0, 1 (42)
The set of differential equations ( 35)-( 37) is an eigenvalue problem where s is the eigenvalue and X 11 = (F 11 , G 11 , G11 ) the eigenvector. It can be written

s M • X 11 = L • X 11 . (43)
Actually, Eq. ( 37) can be solved analytically:

G11 (z) = G 11 (z = 1) sinh k (1 + Λ -z) sinh kΛ ; 1 ≤ z ≤ 1 + Λ , (44)
G11 (z) = G 11 (z = 0) sinh k (Λ + z) sinh kΛ ; -Λ ≤ z ≤ 0 , (45)
with k = k 2 + sκ/κ p . Hence, the eigenvalue problem ( 35)-( 37) can be restrained to the fluid domain, i.e. Eqs ( 35) -( 36), with the boundary conditions

DG 11 = ±ξ kG 11 coth kΛ ; z = 0, 1. (46)
The eigenvalue problem ( 35)-( 36) with the boundary conditions ( 46) is solved using a Chebyshev collocation method. The functions, F 11 and G 11 are expanded in series of Chebyshev polynomial series of order N. The 2 (N + 1) unknowns are determined at the Gauss-Lobatto nodes

z j = 1 2 cos π j N + 1 j = 0, 1, ..., N. (47)
Since k depends on s, an iterative process is implemented. The eigenvector X 11 is normalized such that

G 11 (z = 1/2) = 1. ( 48
)
The marginal stability curve Ra(k) is determined by the condition s = 0. The minimum of the marginality stability curve gives the critical Rayleigh number Ra c and k c respectively. Indeed, when a temperature fluctuation occurs in the liquid close to highly conducting wall, it easily relaxes. Whereas, it can persist and distorts the temperature distribution when the thermal conductivity of the boundary is very low. This temperature distortion can lead to an instability. As a consequence, the critical Rayleigh number is smaller [START_REF] Cerisier | Thermal influence of boundaries on the onset of Rayleigh-Bénard convection[END_REF]. When the convection flow starts, it distorts the originally horizontal isotherms. Since, this deviation from the basic state occurs in the fluid layer as well as in the wall boundary, the wavelength of the convection pattern becomes larger with decreasing ξ [START_REF] Westerburg | Finite amplitude convection in the presence of finitely conducting boundaries[END_REF]. From mathematical point of view, the decrease of Ra c is caused by the weakening of the thermal boundary condition [START_REF] Reynolds | Finite amplitude instabillity of parrallel shear flows[END_REF] for θ as ξ decreases from 10 3 to 10 -3 . For a given thermal conductivities ratio, the critical Rayleigh number for NSBC is greater than that for SFBC. The physical reason is quite intuitive. In the case of NSBC, the friction of the fluid against the wall dissipates more energy, therefore a hogher thermal gradient has to beimposed so that the convection can start.

Additional properties of the critical mode are given by F 11 and G 11 at the critical conditions.

They are displayed in Fig. 2 First, we need to define an inner product between two vectors f and g by

f , g = 1+Λ -Λ f • gdz. ( 49
)
The adjoint mode X ad associated to the critical mode X 11 , verifies the adjoint eigenvalue problem

L + • X ad = 0. ( 50
)
The adjoint operator L + of operator L is defined by

X ad , L • X = L + • X ad , X . (51)

C. Characteristic time of the instability

Near the onset of convection, the growth rate Re(s) of the perturbation may be approximated using Taylor expansion,

s = ǫ τ 0 + O ǫ 2 with ǫ = Ra -Ra c Ra c . ( 52
)
The determination of the characteristic time follows the methodology described by Cross [START_REF] Cross | Derivation of the amplitude equation at the Rayleigh-Bénard instability[END_REF], Plaut [START_REF] Plaut | Modélisation d'instabilités par méthodes nonlinéaires[END_REF]. Explicitly, it is given by

τ -1 0 = -k 2 c Ra c G 11 , F ad G 11 , G ad + G11 , Gad + 1/P r (D 2 -k 2 )F 11 , G ad . ( 53
)
The integrals are evaluated numerically by means of Clenshaw and Curtis method, in terms of the critical conditions. Note that we consider F 11 = G 11 = 0 inside the slabs and G11 = 0 in the fluid domain. Figure 4 shows the variation of τ 0 as a function of ξ. As it can be observed, The characteristic time of instability increases significantly as the thermal conductivities ratio decreases. In other words, the time needed for the system to reach a state of thermal equilibrium increases drastically with decreasing ξ.

IV. WEAKLY NONLINEAR STABILITY ANALYSIS

A. Principles and procedure

A standar weakly nonlinear convection analysis using the amplitude expansion method is adopted as a first approach to investigate nonlinear effects (see Stuart [START_REF] Stuart | On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. part 1. the basic behaviour in plane Poiseuille flow[END_REF], Watson [START_REF] Watson | On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. part 1. the development of a solution for plane Poiseuille flow and for plane Couette flow[END_REF],

Reynolds & Potter [START_REF] Reynolds | Finite amplitude instabillity of parrallel shear flows[END_REF], Sen & Venkateswarlu [START_REF] Sen | On the stability of plane Poiseuille flow to finite-amplitude disturbances, considering the higher-order Landau coefficients[END_REF], Fujimura & Yamada [START_REF] Fujimura | Hexagons and triangles in the Rayleigh-Bénard problem: quinticorder equations on a hexagonal lattice[END_REF], Generalis & Fujimura [START_REF] Generalis | Range of validity of weakly nonlinear theory in the Rayleigh-Bénard problem[END_REF]). At leading order, one writes w(x, y, z, t) = f (x, y, t) F 11 (z) + c.c., (54)

θ(x, y, z, t) = f (x, y, t) G 11 (z) + c.c., (55) θ
(x, y, z, t) = f (x, y, t) G11 (z) + c.c. (56)
The planform function which describes the convection pattern is

f (x, y, t) = N p=1 A p (t) exp (i k p • r) (57)
where |k p | = k c , and A p (t) the amplitude of the perturbation. According to the normalization of the eigenfunctions used in the linear theory, A p (t) represents the amplitude of the thermal perturbation measured at the midplane. Configuration with N = 1 corresponds to rolls and N = 2 to squares. The weakly nonlinear analysis is applied to each of these patterns. The configuration with N = 3 corresponding to hexagons is not considered here.

Further calculations show that this three-dimensional pattern is unstable.

by the fundamental disturbance with wavenumber k = k c , its higher harmonics generated by the nonlinear self-interactions and the modification of the base state due to the interaction with the complex conjugate. As in Stuart [START_REF] Stuart | On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. part 1. the basic behaviour in plane Poiseuille flow[END_REF], Watson [START_REF] Watson | On the non-linear mechanics of wave disturbances in stable and unstable parallel flows. part 1. the development of a solution for plane Poiseuille flow and for plane Couette flow[END_REF], Herbert [START_REF] Herbert | On perturbation methods in nonlinear stability theory[END_REF], the disturbance is expanded in harmonic series and the coefficient of each harmonic is further expanded in an asymptotic series with disturbance amplitude as a small parameter.

In the case of rolls, the velocity and the temperature disturbances are expanded as follows:

w(x, z; t), θ(x, z; t), θ(x, z; t

= m=1 F 0,2m (z), G 0,2m (z), G0,2m (z) |A| 2m + (58) n=1 m=0 F n,n+2m) (z), G n,n+2m) (z), Gn,n+2m) (z) |A| 2m A n E n + c.c. ,
where c.c. means the complex conjugate of its preceding expression, E n = e inkcx ;, n denotes the harmonic index and m indicates the asymptotic order. The time evolution of the amplitude A(t) is given by the Stuart-Landau equation

dA dt = g 0 A + g 1 |A| 2 A + g 2 |A| 4 A + ... (59)
In Eq. ( 59), g 0 = s is the linear growth rate and g 1 is the first Landau constant. The sign of g 1 determines whether the nature of the bifurcation is supercritical (g 1 < 0) or subcritical (g 1 > 0). Substituting expansions (58) and ( 59) into ( 21 It can be shown that g 1 is the sum of contributions of two terms g I 1 arising from the nonlinear inertial terms and g V 1 arising from the nonlinear viscous terms. Since at the lowest order µ = 1 -αΓ, Eqs. ( 13) and [START_REF] Herbert | On perturbation methods in nonlinear stability theory[END_REF], where α = 1-nc 2 λ 2 , then

g 1 = g I 1 + g V 1 with g V 1 -αg N N 1 (60)
In the above sections, the nonlinear stability problem is reduced to a sequence of differential equations. As in the linear problem, they are solved using a spectral collocation method based on Chebyshev polynomials. The differential equations are collocated at Gauss-Lobatto points. The integrals involved in the determination of the first Landau constant are calculated using Clenshaw and Curtis method.

V. RESULTS AND DISCUSSION

A. Bifurcation to rolls

The first Landau constant g 1 as well as the different contributions g I 1 and g N N

1
are determined for different critical sets (Ra c , k c , ξ, P r).

In Fig. 5(a) , we plot g I 1 as a function of ξ. As expected, g I 1 is negative, i.e., the bifurcation is supercritical for a Newtonian fluid. The absolute value of g I 1 decreases with decreasing ξ and |g I 1 | → 0 when ξ → 0. Note that g I 1 is sensitive to change in ξ mainly when 0.1 ≤ ξ ≤ 5. For P r ≥ 1, the analysis of the contribution to g I 1 arising from the different nonlinear interactions shows that g I 1 is dominated by the nonlinear thermal convection terms involving the mofication of the conductive temperature profile, for all the range of ξ considered. The contribution of the nonlinear inertial term is practically negligible.

As it can be observed, g N N 1 is negative, and g

V 1 = -αg N N 1 > 0.
Therefore, shear-thinning effects promote a subcritical bifurcation. In Fig. 6 we plot g 1 as a function of ξ for different values of α. For low shear thinning effects, g 1 < 0 and the bifurcation is supercritical, while for sufficiently high shear-thinning effects, g 1 > 0 and the bifurcation is subcritical. Using Eqs. (60), the critical degree of shear-thinning α c above which the bifurcation changes from supercritical to subcritical is given by

α c = g I 1 g N N 1 . ( 61
)
Figure 7 shows the variation of α c as a function of ξ at P r = 10, for NSBC and SFBC. The more ξ is low, the more the degree of shear-thinning α of the fluid must be high to obtain a subcritical bifurcation. This is due to the reduction of the convection intensity when ξ decreases. This effect is even more pronounced in the case of SFBC. Note that, the largest variation of α c with ξ occurs mainly for 10 -2 ≤ ξ ≤ 10.

Bifurcation to squares

For square patterns, considered as the superposiution of two perpensicular sets of rolls, the amplitude equations are [START_REF] Golubitsky | Symmetries and pattern selection in Rayleigh-Bénard convection[END_REF]:

dA 1 dt = s A 1 + g 1 A 2 1 + λ 1 A 2 2 A 1 , (62)
dA 2 dt = s A 2 + g 1 A 2 2 + λ 1 A 2 1 A 2 . (63)
As for g 1 , the coefficient of cross saturation λ 1 is obtained by invoking the solvability condition and can be written as: The numerical results indicate that λ I 1 and λ N N 1 are negative and their variation with ξ is similar to that of g I 1 and g N N 1 respectively. The critical value of α above which, the bifurcation becomes subcritical is given by

λ 1 = λ I 1 -αλ N N 1 . (64) 10 -3 10 -2 10
α c = g I 1 + λ I 1 (g N N 1 + λ N N 1 ) (65)
The variation of α c as a function of ξ, at P r = 10, is depicted in Fig. 8 for NSBC and SFBC. As in the case of rolls, α c increases with decreasing ξ.

The influence of Prandtl number on α c is shown in Figure 9 for two limit values of ξ:

10 -3 and 10 3 . Note that for P r ≥ 10, α c no longer depends on P r. However, when P r < 1, α c increases strongly as P r decreass.

B. Pattern selection

In this section, we investigate the pattern selection on a square lattice. The calculation proceeds in two stages. First, the possible steady states solutions of the amplitude equations are determined. Then, their linearized stability is determined by computing the eigenvalues of the linearized system around each solution. A solution is stable if all its eigenvalues are negative.

The system of amplitude equations for a square lattice are given by Eqs (62) and (63). The coefficients s, g 1 and λ 1 in these equations depend on the rheological parameters and on the reduced Rayleigh number ǫ = (Ra -Ra c ) /Ra c . The stationary solutions are obtained by

setting f i (A 1 , A 2 ) = 0
, where f i is the right hand side of the amplitude equations. Their stability is determined by the sign of the eigenvalues χ i of the Jacobian matrix J ij = ∂f i ∂A j evaluated at the steady states. In the following, the stability of the stationary solutions is examined in details.

(i) Conduction state, A 1 = A 2 = 0. The eigenvalues associated to this state are χ 1 = χ 2 = s.

The conduction state is stable if ǫ < 0 and undergoes a stationary bifurcation at ǫ = 0.

(ii) Steady convection with rolls parallel to e x or e y , A 1 = -s/g 1 , A 2 = 0 or -s/g 1 . The eigenvalues associated to this state are: χ 1 = -2s and χ 2 = s g 1 -λ 1 g 1 .

A 1 = 0, A 2 = 10 -1 10 0 0.8 1 1.2 ξ λ 1 /g 1 (1) (4)
(iii) Steady convection with square patterns, A 1 = A 2 = -s/(g 1 + λ 1 ). The eigenvalues associated to this steady state are χ 1 = -2s and χ 2 = 2 s (λ 1 -g 1 ) λ 1 + g 1 .

In the supercritical regime, i.e. s > 0 and χ 1 < 0. The sign of χ 2 depends on the ratio λ 1 /g 1 . α = 0, ξ c = 1. This result is in agreement with that given by Jenkins and Proctor [START_REF] Jenkins | The transition from roll to square-cell solutions in Rayleigh-Bénard convection[END_REF].

It is represented in

With increasing shear thinning effects, the interaction between the two sets of rolls of a square-cell increases, reducing by this way, ξ c . order become large and should be taken into account. A weakly nonlinear expansion was then carried out up to fifth-order in amplitude. Figure 12 shows the evolution of ξ c versus the reduced Rayleigh number, ǫ, for different values of the constant time of the fluid λ. The shear-thinning index is fixed to n = 0.5. The intensity of convection increases with increasing ǫ. The interaction between the two sets of rolls, via nonlinear inertial and nonlinear viscous terms become stronger leading to a diminution of ξ c .

VI. SOLUTIONS AT HIGHER ORDER

In the present section, information on the heat transfer, the flow structure and viscosity field in rolls and squares are provided. The influence of shear-thinning will be emphasized.

A. Heat transfer

The heat transfer through the horizontal fluid layer is described by the Nusselt number, Nu, the ratio of the total heat flux to the purely conductive heat flux in the absence of fluid flow. It can be calculated either at the lower or upper plate. At the lower plate, we have

Nu = 1 - ∂ θ ∂z z=0 = 1 - N p=1 M m=1 A 2m p (DG 02m ) z=0 -S A 2 1 A 2 2 (DG 04 ) z=0 (66)
where the overbar denotes the horizontal average over one wavelength, N = 1 corresponds to rolls and N = 2 to squares, M = 1, when the series ( 59) is truncated at the third order, and M = 2, when ( 59) is truncated at the fifth order, S = 1 for squares and S = 0 for rolls.

The term DG 04 arises from the interaction between modes with different eigenvectors. The unperturbed solution, Nu = 1, corresponds to the hydrostatic solution. The second term of Nu refers to the convective transfer.

Figure 13 shows for a two-dimensional roll solution computed at the fifth-order, the evolution of Nu -1 as a function of ǫ for different values of ξ. The Nusselt number decreases with decreasing ξ. When ξ → 0, i.e. for poorly conducting plates, Nu → 1. The convection in the fluid layer will not contribute to the overall heat transfer. For ξ = 1000, to represent the case of perfectly conducting walls, our results are in good agreement with the numerical solution of ( 21) and ( 22) obtained using the spectral code of [START_REF] Plaut | Low-Prandtl-number convection in a rotating cylindrical annulus[END_REF], at least up to ǫ = 0. & Veronis [START_REF] Malkus | Finite amplitude cellular convection[END_REF])

B. Viscosity field

The viscosity distribution for rolls and square cells, computed at the fifth order in amplitude, is shown in Fig. 15 For ξ > 10, α c is practically the same as for perfectly conducting plates. The numerical results indicate also that α c is independent of P r when P r > 10. The stability of rolls and squares is then investigated as a function of ξ and the rheological parameters. In the Newtonian case, squares are stable when ξ < 1, in agreement with [START_REF] Jenkins | The transition from roll to square-cell solutions in Rayleigh-Bénard convection[END_REF]. In the case of shear-thinning fluids, an additional nonlinear coupling between modes is introduced by the rheological law. This leads to a decrease of the crtical value of ξ below which squares are stable.

By considering the amplitude expansion at the fifth order, the range of validity of the weakly nonlinear analysis is extended and the domain of stability of the square pattern in the (ξ, ǫ) plane is determined. These results are consistent with the maximum heat-transfer principle.

In the weakly nonlinear approach, it is assumed implicitly that the dynamics is dominated by the fundamental mode. We intend to analyse the stability of the convective patterns, in

Figure ? ?

 ? Figure ?? displays the variation of Ra c and k c as a function of the ratio ξ of the thermal conductivity of the boundary to that of the fluid in the case of NSBC and SFBC. The dimensionless thickness of the plates is fixed at Λ = 1. These results are in very good quantitative agreement with those obtained by Sparrow et al.[START_REF] Sparrow | Thermal instability in a horizontal fluid layer: effect of boundary conditions and a non-linear temperature profile[END_REF], Proctor[START_REF] Proctor | Planform selection by finite amplitude thermal convection between poorly conducting slabs[END_REF], Jenkins and Proctor[START_REF] Jenkins | The transition from roll to square-cell solutions in Rayleigh-Bénard convection[END_REF], Carriere et al.[START_REF] Carriere | Wavelength selection in Rayleigh-Bénard convection between horizontal boundaries of finite conductivity[END_REF] and Cerisier et al.[START_REF] Cerisier | Thermal influence of boundaries on the onset of Rayleigh-Bénard convection[END_REF]. The wavelength of the convection becomes larger with decreasing ξ and the critical Rayleigh number Ra c is also reduced.

10

 10

FIG. 1 .

 1 FIG. 1. Critical Rayleigh number (a) and critical wavenumber (b) as function of ξ in the case of NSBC (1) and SFBC (2).

  for NSBC and SFBC and different values of ξ. As indicated above, when ξ decreases, the temperature perturbation persists at the wall (Figs.1c 1d, curves 3, 2 and 1), reducing Ra c and the velocity perturbation is damped (1a and 1b, curves 3, 2 and 1). The coupling betweeb the temperature perturbation θ in the fluid and that in the horizontal plates θ is illustated by Fig.3(a), where contours of the temperature perturbation at ξ = 1 are represented. For comparison, the case of large value of ξ, where the temperature vanshes in the liquid vanishes at the boundaries is also represented (3 a).

FIG. 2 .

 2 FIG. 2. Eigenfunctions at critical conditions and different values of ξ: (1) ξ = 10 -3 ; (2) ξ = 10 -1 ; (3) ξ = 10 0 and (4) ξ = 10 3 in the case of NSBC (a-c) and SFBC (b-d).

FIG. 3 .

 3 FIG. 3. Temperature perturbation contours for two values of the thermal conductivity ratio ξ: (a) ξ = 1 and (b) ξ = 1000. In this latter case, the temperature perturbation in the liquid relaxes at the boundaries.

FIG. 4 .

 4 FIG. 4. Characteristic time of instability vs thermal conductivities ratio ξ at Pr = 10. (1) SFBC and (2) NSBC.

  )-(23) yields after some algebra the differential equation for any F n,2m+n , G n,2m+n and Gn,2m+n , which are solved sequentially beginning from n = 1 and m = 0. The problem with harmonic index n = 1 and amplitude order m = 0 is the linear stability problem (35)-(37). The problem n = 0, m = 1 is O(A 2 ) correction of the conductive temperature profile due to nonlinear interactions. The problem n = 2, m = 0 is the first harmonic of the fundamental mode which manifests at order O(A 2 ). The problem n = 1, m = 1 is O(A 3 ) correction to the fundamental mode. It contains nonhomogeneous terms due to nonlinear interations. The application of the Fredholm alternative allows the determination of the first Landau constant (the first correction to the linear growth rate) which appears in the time derivative of w, θ and θ.

FIG. 5 .

 5 FIG. 5. Contribution of the nonlinear inertial terms (a) and nonlinear viscous terms (b) to the first cubic Landau constant, as a function of the ratio ξ at P r = 10. Curves (1) and (2) correspond to NSBC and SFBC respectively.

FIG. 6 .

 6 FIG. 6. Cubic Landau constant as a function of ξ in the case of NSBC with P r = 10 and different values of α: (1) α = 0, i.e. Newtonian case ; (2) α = 1 × 10 -5 ; (3) α = 5 × 10 -5 ; (4) α = 10 -4 ; (5) α = 2.5 × 10 -4 ; (6) α = 3.5 × 10 -4 ; (7) α = 5 × 10 -4 .

FIG. 7 .

 7 FIG. 7. For rolls, critical value of the degree of shear-thinning α c , vs the thermal-conductivities ratio ξ at P r = 10 for (a) NSBC and (b) SFBC.

FIG. 8 .

 8 FIG. 8. For squares, critical value of the degree of shear-thinning α c , vs the thermal-conductivities ratio ξ at P r = 10 for (a) NSBC and (b) SFBC.

FIG. 9 .

 9 FIG. 9. Critical value of the degree of shear-thinning α c = dµ dΓ Γ=0 c as a function of Prandtl number. Case of square patterns. (1) NSBC with ξ = 10 3 ; (2) SFBC with ξ = 10 3 ; (3) NSBC with ξ = 10 -3 ; (4) SFBC with ξ = 10 -3 .

( 5 )FIG. 10 .

 510 FIG. 10. Ratio of the cross saturation coefficient λ 1 to the saturation coefficient g 1 versus the ratio of the thermal conductivities ξ for different values of the shear-thinning degree α in the case of NSBC. (1) α = 0; (2) α = 2.5 × 10 -5 ; (3) α = 5.625 × 10 -5 ; (4) α = 10 -4 ; (5) α = 1.5625 × 10 -4 .

  Fig. 10 as a function of ξ for different values of α. When |λ 1 | > |g 1 |, χ 2 is positive and the squares are unstable. According to [10], the interaction between the two sets of rolls is too strong and one of the two sets of rolls nonlinearly damps out the other. When |λ 1 | < |g 1 |, χ 2 < 0 and the squares are stable. The critical value of ξ at which the planform of convection changes from square-cell solution (ξ < ξ c ) to two-dimensional roll solution (ξ > ξ c ) is given as a function of α in Fig. 11. In the Newtonian case, i.e.

Figure 11 FIG. 11 .FIG. 12 .

 111112 Figure 11 is obtained by truncating the series (59) to the first Landau constant, i.e. at cubic order in A. For a significant deviation from the critical conditions, terms of higher

FIG. 13 .

 13 FIG. 13. Nusselt number as a function of the reduced Rayleigh number ǫ in the case of NSBC at P r = 10 and for three values of ξ. The Continuous lines correspond to a Newtonian fluid and the dashed lines to a shear-thinning fluid with n = 0.5, λ = 0.02; ( ) numerical solution of (21) and (22) obtained using the spectral code of Plaut & Busse[START_REF] Plaut | Low-Prandtl-number convection in a rotating cylindrical annulus[END_REF] in the case of a Newtonian fluid with perfectly conducting walls.

FIG. 14 .

 14 FIG. 14. Nusselt number as a function of the reduced Rayleigh number ǫ in the case of NSBC with P r = 10 and for two values of ξ. (a) Newtonian fluid, (b) shear-thinning fluid with n = 0.5 and λ = 0.02. The Continuous lines correspond to rolls and the dashed lines to squares.

FIG. 15 .

 15 FIG. 15. Flow structure and viscosity distribution over a roll (a, ξ = 1000) and a square (b, ξ = 0.1) at ǫ = 0.2, with NSBC, α = 0.2α c and P r = 10. For the square, the horizontal velocity field is shown at z = 0.8. (c) Vertical velocity profile in the square cell centre for NSBC and SFBC.

.

  figure), the viscosity is weakly reduced. The viscosity is minimal at locations (dark region in the figure) where the shear rates γxz and γyz are maximal.

Acknowledgements

The authors gratefully acknowledge the financial support of the French National Research Agency, under the grant ANR-10-BLAN-925-01.


L Pizzella
email: laure.pizzella@mines-paristech.fr

C Lajaunie

G Courrioux

S Lopez

X Freulon

D Renard

B Bourgine

Integrating field interpretations to geological modeling with the potential field method

Geological structures, such as complex folds, fluvial channels or ore-veins, usually present principal directions of continuity that can be defined locally but vary on a broader scale. Taking into account both these local anisotropies and their variability is key to achieve sound geological 3D modeling. We make the hypothesis here that anisotropy variations can be sampled on the field and incorporated into the model at an early stage of model building. This work develops two strategies to incorporate such observations to implicit modeling using the potential field method: i) adding second-order derivative information to this field; ii) describing the field as a stochastic convolution in which will be given an interpolated anisotropic field.

Introduction

Over the last decades, new geological modeling methods have proven to be an efficient approach to quickly build complex structural models. They are based on the description of geological interfaces as implicit surfaces and their boolean combination, using Constructive Solid Geometry rules with a geological meaning. Among them, the so-called potential field method [START_REF] Lajaunie | Foliation fields and 3D cartography in geology: principles of a method based on potential interpolation[END_REF][START_REF] Calcagno | Geological modelling from field data and geological knowledge, Part I -Modelling method coupling 3D potentialfield interpolation and geological rules[END_REF] is based on the co-kriging of field data in order to construct a scalar field Z. This way, a surface is designed as a particular isovalue of the field in 3D space. Those data can derive from drill-hole data, digital geological maps, structural data, interpreted cross-sections and consist in:

• contact points, sampling a same geological interface k: Z(u j ) -Z(u j ) = 0 for all pair (j, j ) of any set J k

(1)

• orientation data such as dip directions or tangents, expressed using the gradient of the scalar field:

∂Z(u i ) ∂x = G x i (2)
The interpolation of the scalar field is done thanks to a co-kriging:

[Z(u) -Z(u 0 )] * = k jj ∈k λ k,j (Z(u j ) -Z(u j )) + i dim t β i,t G t i (3)
Though such an approach has provided very satisfactory results, available data is frequently not sufficient to describe complex structures, especially in cases of locally varying anisotropies [START_REF] Martin | Non-stationary Geostatistical Modeling: A Case Study[END_REF]. When facing important uncertainties, experienced users have to introduce a priori geological knowledge into the model in the form of additional control points. This workaround is usually time consuming, error prone and user dependent. It may also lead to overconstrained models which no longer leave room to uncertainty analysis and/or quantification. This work will present ways to incorporate geological observations of anisotropy at an earlier stage of modeling. The first one is based on second-order derivatives data of the scalar field and the second one on a description of the potential as a stochastic convolution of a Gaussian kernel.

Contribution of second-order derivative information

Just as the gradient is an information of the first-order derivative of the potential field, the second-order derivative can be integrated into the co-kriging equations. They represent the curvature of the field at a point u i in the specific direction d:

∂ 2 Z(u i ) ∂ d 2 = lim →0 Z(u i + d) + Z(u i -d) -2 Z(u i ) 2 d 2 (4)
It could seem hard to sample in practice but actually special cases of null second-order derivatives are already measured, by fold axis data for example. However, whereas introducing gradients of the potential field requires a twice differentiable covariance function K, considering covariance between second-order derivatives requires the function to be differentiable four times:

< ∇ 2 Z(u), ∇ 2 Z(v) > = ∇ 4 u,u,v,v K(u, v) = ∇ 4 u,u,v,v C(r) where r = (u -v) t A(u -v) (5)
This way, the cubic covariance traditionally used for the potential method no longer applies. One covariance function which fulfills the requirements and then considered in this work is a radial one, piecewise polynomial, positive definite and compactly supported [START_REF] Wendland | Piecewise polynomial, postitive definite and compactly supported radial functions of minimal degree[END_REF]: This approach can integrate fold axis, which are real field data, traditionally described but unused in practice. The contribution of these second-order derivatives values is illustrated on a 2D synthetic example, shown in Figure (1).

C(r) = 1 3 (1 -r) 6 + (35r 2 + 18r + 3) (6)
Another intended approach to integrate anisotropy is to model the potential as a convolution product between a stochastic white noise W and a deterministic kernel Φ, slowly varying with the position of u (Higdon, 1998;[START_REF] Fouedjio | A generalized convolution model and estimation for non-stationary random functions[END_REF]:

Z(u) = R Φ Bu (u -s) W (ds) (7)
This kernel Φ requires a complete definition of the anisotropy since it defines at each point u, a local anisotropy B u , according to 4 (respectively 7) functions in 2D (respectively 3D):

• θ(u) the main rotation angle of the anisotropy direction in 2D, α(u), β(u) and γ(u) the 3 rotation angles of the main anisotropy direction in 3D

• a 1 (u) and a 2 (u), the scale factors of anisotropy in respectively main and orthogonal direction of anisotropy, a 1 (u), a 2 (u) and a 3 (u) in 3D

• σ(u), the standard deviation A Gaussian kernel has been chosen since its derivative, which are required for derivations of Z, are explicit. Figure (2) presents a 2D example of stochastic convolution use on the same case than Figure (1). The function of the angle of local anisotropy θ(u) is interpolated thanks to a complex kriging [START_REF] Lajaunie | Sur le krigeage des fonctions complexes[END_REF] of some directional data.

Conclusion and perspectives

This work presents two different approaches in order to take into account issue of non-stationarity of some structural implicit modeling cases. Both can rely on real field data as well as interpretative descriptions. Whereas the second-order derivative method may appear as a local description of anisotropy, the stochastic convolution method is a more global approach but in return requires an explicit function of the anisotropy angle. Thus one key of this approach is the angle interpolation, performed here with complex kriging in 2D, but which could eventually be achieved in 3D with quaternion interpolation. These methods are not to be opposed, since each can benefit from the other: it is theoretically feasible -even though it makes the calculations harder -to take into account the second-order derivative data in the convolution method. In future work, these improvements have to be challenged on real 3D datasets which present such variable anisotropy.

Figure 1 :

 1 Figure 1: Impact of the second-order derivative data into the potential field method. Red points (6) are contact points (increments) of a same iso-potential. Black arrows (6) are gradient values of the field. The actual result (a) can not reproduce the continuity of the structure. Addition of second-order derivative values, constrained to zero, in this direction ((b), blue arrows) allows the method to follow the interpreted anisotropy.

Figure 2 :

 2 Figure 2: Stochastic convolution method results