Emilien Sauvage
email: emilien.sauvage@hmg.inpg.fr

Laetitia Jacoutot

Yves Fautrelle

Annie Gagnoud

Laure Blumenfeld

Eric Favre

François Daviaud

Comparison between numerical and experimental results on thermoconvective instabilities of a high-Prandtl-number liquid

Keywords: number͑s͒: 47, 20, Ϫk, 47, 35, Ϫi

The flow structuration of silicon oil ͑Prandtl number of 10.3͒ in a open cylindrical pool heated from the center of the surface is investigated numerically. Our purpose is to perform the numerical simulation of experimental results obtained by Favre et al. ͓Phys. Fluids 9, 1473 ͑1997͔͒ who observed transitions between steady and axisymmetric flows at sufficiently low values of the Marangoni number ͑Ma͒ and various types of instability depending on the height of the fluid. The hydrothermal wave regime has been obtained at critical values of Ma which depend on the Bond number and on the aspect ratio. The numerical results are in good agreement with the experimental ones.

I. INTRODUCTION

When a liquid free surface is submitted to a horizontal thermal gradient, there appears a surface stress which drives a flow along the interface. The latter phenomenon called the thermocapillary effect, also called the Marangoni effect, comes from the existence of variations of the surface tension along the interface correlated to the temperature variations ⌬T. Generally, thermocapillary effects are coupled with natural convection. Since buoyancy is a body force, its importance depends on the volume of the pool in comparison with its free surface.

The thermocapillary effect is found in various industrial processes. For example, in welding devices by means of electron/laser beam or arc, the molten pool is heated by its free surface, and there exist very strong temperature gradients along the surface ͓1-3͔. The Marangoni effect is usually dominant, and significant fluid flows are observed. Similar situations are encountered in evaporating devices which use concentrated heat sources ͑laser or electron beam͒ for the coating of substrates by vapor deposition.

In crystal growth, thermocapillary effects are also present, although less intense ͓4͔. Nevertheless, Marangoni convection is responsible for various defects, e.g., striations of the solid crystal, which result from the flow instabilities of the thermocapillary type.

The aim of the present paper is to analyze, thanks to numerical calculations, the various types of flow configurations obtained experimentally by Favre ͓2,3͔. First the different types of flow and instabilities put forth in the literature will be reviewed hereafter in the present section. Then in Sec. II the formulation of the problem is provided. In Sec. III the various types of computed flow patterns and instabilities are compared to the experimental data.

Basic flow. Contrary to the Bénard type flow ͑vertical thermal gradient only͒ for which a threshold of temperature exists, in this configuration as soon as ⌬T 0, the fluid moves. For example, if we consider the simplest case corresponding to a fluid cavity submitted to a horizontal temperature difference ⌬T applied on each lateral side, the flow near the equilibrium state, i.e., small ⌬T, is bidimensional and stationary. This flow is named "basic flow." Hydrodynamic theory for the thermocapillary effect is quite old since the basic flow in a rectangular cavity with no gravity was solved by Birikh ͓5͔ in 1962 and the first linear stability study was made in 1983 by Smith and Davis ͓6͔. Laure et al. ͓7͔ and Mercier and Normand ͓8͔ have extended this solution to a flow with thermogravity effect and other boundary conditions. Finally, Garnier ͓9͔ proposed an analytical solution to the basic flow in a cylindrical cavity based on a method introduced by Laure et al. ͓7͔. This solution will be extended to the specific case studied experimentally by Favre ͓2͔ experiments in which the bottom boundary is cooled.

Instabilities. When the imposed temperature difference ⌬T crosses a certain threshold, the basic flow is destabilized. Instabilities have various forms and characteristics depending on many parameters of the problem such as the Prandtl number, Bond number, aspect ratio, Biot number, and the type of boundary conditions. The most studied case in the literature is a rectangular cavity heated and cooled on both sides by a vertical boundary, the free surface and the bottom being insulating or conducting. In this configuration the thermal gradient is horizontal, which simplifies the analytical solution. Annular geometry is more and more studied because of two reasons. First, in this particular case, the absence of walls in the orthoradial direction ͑i.e., the direction perpendicular to the thermal gradient͒ avoids experimentally wave reflection on the side walls. Second, this annular ge-ometry is closer to the industrial process of the Czochralski crystal growth method of electron beam and laser surface heating. Very few studies ͓10,11͔ have considered a heating configuration localized on the surface, as used in the Favre ͓2͔ study. Moreover, the condition of a cooled bottom boundary is also rarely used ͓12͔. The strong variation of the stability criteria with the numerous problem parameters makes the comparison between results of different studies very hard to achieve. However, the instabilities can be sorted out in three groups.

The first type of instability in a similar configuration is the function of the aspect ratio of the fluid layer. For a small Bond number, this instability is composed of corotative stationary or time dependent rolls with a two-dimensional longitudinal axis. Laure et al. ͓7͔ give critical values associated to their apparition in a long rectangular cavity with a low Prandtl number fluid. Mercier and Normand ͓13͔ have shown that rolls appear near the cold wall for fluid with a low Prandtl number ͑PrϽ 0.01͒, near the hot wall for fluid with a high Prandtl number ͑PrϾ 4͒, and near both sides for fluid with an intermediate Prandtl number. These results have been checked in a great number of two-dimensional numerical publications ͓14-18͔. The influence upon the flow structure of parameters as the aspect ratio, the Bond number, and several boundary condition types has been also studied. Some numerical studies ͓14-16͔ report the presence of temporal oscillations due to the destabilization of the first corotative roll. Experimentally this instability is well observed for various configuration ͓5,12,15,19͔.

The second type of instability consists of stationary rolls with an axis collinear to the thermal gradient direction. This instability has been reported experimentally by Daviaud and Vince ͓20͔ in a rectangular cavity and in an annular cavity by Favre ͓3͔ for a relatively high height of fluid ͑h Ͼ 3 mm͒. Li et al. ͓11͔ found numerically longitudinal rolls in an annular cavity while Mercier and Normand ͓8͔, from an analytical point of view, have predicted this instability for a fluid of Prandtl number Pr= 7 in a cavity with conducting thermal boundary conditions at the bottom and introducing a Biot number at the surface. This kind of instability has not been studied widely.

The third and most interesting instability are the hydrothermal waves. Smith and Davis ͓6͔ have studied in 1983, the linear stability of a return flow without gravity and discovered this new kind of instability. They appear as a propagative ondulatory instability of the basic flow with specific frequency and wave number. The propagation angle of the waves with the thermal gradient direction is a function of the Prandtl number ͓6͔. Numerous experiences observed these hydrothermal waves in rectangular ͓19͔ and annular ͓12͔ geometry. A summary is presented by Burguete et al. ͓21͔. A stability diagram of the flow is achieved plotting the critical values of imposed thermal difference ͑or Marangoni number͒ for which instability occurs against the height of the fluid in the cavity ͑or Bond number͒. The critical Marangoni value increases with the height of fluid. If the height of fluid crosses a certain threshold the instability is not hydrothermal waves any longer but an oscillating multicellular flow. In a majority of studies, the hydrothermal waves are the second instability, the first being the corotative instability. A similar stability diagram of the flow has been drawn by Favre ͓2͔ which is shown in Fig. 1. Garnier ͓9͔ puts forth the assumption that the hydrothermal waves can appear only when Bd Ͼ 1, i.e., when thermogravitational forces are negligible, a suggestion which is in agreement with the theory of Smith and Davis ͓6͔ where only the thermocapillary force was taken into account.

From a numerical point of view, many calculations, especially in two-dimensional geometries, has been achieved ͑see, for example, Ben Hadid and Roux ͓17͔͒. The increase of computation capacity leads to the achievement of simulations of the three-dimensional hydrothermal waves. Hoyas and co-workers ͓22,23͔ performed a series of threedimensional calculations in an annular geometry with a numerical method based upon the resolution of the set of stability equations with a decomposition of the fields in Chebyshev polynoms. They found, for a fluid with an infinite Prandtl number, hydrothermal waves propagating near the surface and studied the threshold of appearance depending on the aspect ratio and thermal boundary. Li and co-workers ͓11,24-26͔ simulate hydrothermal waves in an annular geometry heated either by the outer wall or the surface, solving directly the Navier-Stokes equations with a volume-of-fluid method ͑Pr 1͒. Specifically, they were able to determine both the threshold of apparition of the waves and to visualize them correctly. In 2006, Shi and Imaishi ͓27͔ simulated, with the same methods, a fluid of Prandtl number Pr= 6.7. Visualization of the hydrothermal waves computed by Shi and Imaishi ͓27͔ is similar to that obtained experimentally by Garnier ͓9͔.

Experimental results of Favre ͓2͔. The studied geometry is original because the fluid is heated on the surface in the center of the cavity and cooled both by the bottom and side walls. cillating rays ͑ORs͒. The structures named "rays" correspond to the second type of instabilities described below, i.e., stationary rolls with an axis collinear to the thermal gradient direction. The "petal" structure has been described in the literature as the superposition of two counter-rotative wave trains but Favre ͓2͔ was unable to describe it precisely.

II. PROBLEM FORMULATION

A. Physical and mathematical models

The problem consists of the analysis of a flow in a cylindrical pool with a free surface subjected to a thermal gradient ͑cf. Fig. 3͒. The fluid layer of height h is heated at the surface in the center of the cylindrical pool with a heater of radius R 1 = 15 mm and cooled at the bottom and side of the pool of external radius R 2 = 50 mm. The radial length of the pool is L = R 2 -R 1 , the heater is at the temperature T H , while the cooled surfaces are at the temperature T c . The typical thermal gradient in the pool is ⌬T / L = ͑T H -T c ͒ / L Ն 0. We introduce a Biot number Bi to characterize convective heat transfer on the free surface. The driving forces of the flow are ͑i͒ the buoyancy force which tends to drive the hot fluid upward, and ͑ii͒ the thermocapillary force which draws the surface fluid from the hot central zone toward the cold boundary. The thermocapillary force is generated by the variation of the surface tension with the temperature which is well represented by the linear law ͑T͒ = 0 -␥͑T -T 0 ͒ where ␥ ץ-= / ץT is for common fluids always positive and constant and 0 is the value of the surface tension at the temperature T 0 . The physical properties of the silicon oil are summarized in Table I.

Governing equations. The incompressible governing equations are expressed in a cylindrical reference frame and take into account the buoyancy effect with the Boussinesq approximation,

0 ͫ ץv ជ ץt + ͑v ជ • ٌ ជ ͒v ជ ͬ = -ٌ ជ P + ᭝ v ជ + 0 ␤͑T -T 0 ͒g ជ, ٌ ជ • v ជ = 0, ץT ץt + v ជ • ٌ ជ T = ᭝ T, ͑1͒
where v ជ is the velocity vector such as v ជ͑r , , z͒ = ue r ជ + ve ជ

+ we z ជ , ␤ = 1 0 ץ
ץT is the thermal expansion coefficient supposed to be constant, and 0 is the fluid density at the temperature T 0 .

The boundary conditions. At the bottom z = 0, we use the no-slip condition v ជ͉ z=0 =0 ជ , while the temperature is imposed T͉ z=0 = T c . For z = h, the free surface is assumed to be flat and nondeformable, and a no-flux condition is applied, v ជ • n ជ =0. Furthermore, the viscous surface stress is equal to the capillary stress induced by the thermal gradient:

͉ ץu ץz ͉ z=h = ␥͉ ץT ץr ͉ z=h . The thermal condition is ͉ ץT ץr ͉ z=h +Bi͑T -T ϱ ͒ =0
, where Bi is the Biot number and T ϱ is the temperature far from the surface being equal to T cond . We introduce the nondimensional parameter B such that B =Bi/ ͑1+Bi͒.

This system is characterized by six dimensionless parameters: ͑i͒ Prandlt number: Pr= / ; ͑ii͒ Rayleigh number Ra = g␤⌬Th 4 / L; ͑iii͒ Marangoni number Ma= ␥⌬Th 2 / 0 L; ͑iv͒ aspect ratios = h / L and ͑v͒ = L / R 1 ; ͑vi͒ Biot number Bi= h N h / where h N is the convective heat transfer coefficient at the free surface. Instead of the Rayleigh number we use the dynamic Bond number Bd= Ra/ Ma. The dynamic Bond number can be written as the square of the ratio between the height of fluid h and the capillary length h e 2 = ␥ / 0 g.

B. Numerical methods

The commercial software FLUENT is used for the simulations. The equations are discretized by the finite volume method. The third-order MUSCL scheme ͑monotone upstream-centered schemes for conservation laws͒ is used for the convective terms while a second-order centered Temperature coefficient of surface tension

␥ =-8ϫ 10 -5 N m -1 K -1 Thermal expansion coefficient ␤ = 1.34ϫ 10 -3 K -1
Prandtl number Pr= 10.3 scheme is applied to the diffusion terms. The pressure terms are interpolated with the PRESTO! scheme. The SIMPLEC algorithm is used to handle the pressure coupling. The high orders of resolution are useful to simulate correctly the strong coupling between the flow and the temperature field due to the high value of the Prandtl number ͑Pr= 10.3͒. Convergence at each time step is considered when residual error of the continuity equation and velocity component became less than 10 -5 and less than 10 -7 for the temperature. Numerical simulations were performed on an Intel Bi-pro Xeon 3.80 Ghz, 3.89 Go RAM.

Nonuniform staggered grids are used in this study. The mesh is finer near the free surface and near the cold vertical boundary because of the strong velocity and thermal gradient due to the Marangoni effect and the thermal boundary layer. To assess the effects of numerical diffusion we use the mesh Peclet number, Pe m , which is the ratio of a convection time to a thermal diffusion time at the scale of the mesh, Pe m = u m l m / where u m is an average velocity and l m is a mesh characteristic length. A high value of the Peclet number means an important numerical diffusion and gives less accurate results. Several grids were tested. The convergence tests for two-dimensional grids are presented in Table II for the cases h = 14.9 mm and 2.9 K, while Table III shows the 3Dgrid convergence test for h = 2 mm. The criterion are ͑i͒ the value of the heat flow evacuated by the heater, ͑ii͒ the radial velocity value at the free surface for three different radii ͑25, 30, and 38 mm͒, ͑iii͒ the maximal and global average value of the Peclet number, and ͑iv͒ mesh percentages of the field having a Peclet number, lower than 1 and higher than 10.

It must be noticed that the values of the radial velocity on the surface are very sensitive to the mesh quality, but the global flow structuration does not really change.

The shadowgraphic method used experimentally by Favre consists of averaging the horizontal thermal gradient along the height of fluid. From a numerical point of view, similar visualizations are obtained by deleting the visual effect of the main radial thermal gradient to emphasize azimuthal variations of the thermal gradient. A first method due to Li et al. ͓24͔ consists in computing a fluctuating temperature field TЈ defined as

TЈ͑r,z͒ = T͑r,,z͒ - 1 2 ͵ 0 2
T͑r,,z͒d. ͑2͒

The second method consists in post-treating directly the radial or tangential temperature gradient field.

III. RESULTS AND DISCUSSIONS

A. Basic flow

The analytical solution of the thermoconvective basic flow for the cylindrical geometry within the assumption of small aspect ratio is presented in this section. This solution, developed by Garnier ͓9͔, generalizes those obtained for rectangular geometry. We extend it to the specific experimen- The governing equations for the fluid layer are expressed in a cylindrical reference frame, but to approach the rectangular solution, a change of variable is done, r = R 1 + x. The basic flow being axisymmetric, it implies ץ ץ ϵ 0, v ϵ 0. Henceforth, in all the equations the variables will appear in nondimensional form notified by a star superscript. The scalings are shown in Table IV. The method is based on the expansion of all unknown variables in terms of the aspect ratio considered as a small parameter = h / L, as follows:

u * ͑x * ,z * ͒ = ͚ i=0 n i u i * ͑x * ,z * ͒ + O͑ n ͒.
The above expansion is introduced in the governing equations. The zeroth-order system can be solved analytically. The temperature solution exhibits a linear z dependence because cooled bottom boundary

T 0 * ͑x * , z * ͒ = B T cond * ͑x * , z * ͒,
where the nondimensional conducting profile of temperature T cond * ͑x * ͒, is defined as

T cond * ͑x * ͒ = 1 - ln͑1 + x * ͒ ln͑1 + ͒ , = L / R 1
characterizing the curvature of the cavity. We introduce ⌽͑x * ͒ as the horizontal temperature gradient of the conducting temperature profile,

⌽͑x * ͒ = dT cond * dx * = - 1 ln͑1 + ͒ 1 + x * .
The zeroth-order vertical velocity v 0 * vanishes, while the horizontal component is

u 0 * ͑x * ,z * ͒ = B ⌽͑x * ͒ Pr ͩ Ra 240 ͑10z *4 -27z *2 + 14z * ͒ - Ma 4 ͑3z *2 -2z * ͒ͪ . ͑3͒
In the case of a rectangular geometry, terms of order equal or greater than zero are all vanishing, so an exact solution for the basic state can be found. The curvature of the cavity implies a x dependence of all variables which induces nonzero first-, second-, and higher-order terms. The system of first-order equations, thanks to the zeroth-order solution, can be solved in the same way. The first-order component of temperature is

T 1 * ͑x * ,z * ͒ = B 2 ⌽ 2 ͑x * ͒ 240 ͩ Ra 420 ͑100z *7 -567z *5 + 490z *4 + 175z * -198 B z * ͒ -Ma͑9z *5 -10z *4 -5z * + 6 B z * ͒ͪ . ͑4͒
The expression of the first-order horizontal velocity component u 1 * is determined by the equation

ץ 3 u 1 * ץz *3 = Ra Pr ץT 1 * ץx * +
ץ ץz ͩu 0 * ץu 0 * ץx * ͪ.

Results for the basic flow. The theoretical results developed in the previous section are valid for a small aspect ratio, typically Յ0.05. In the studied configuration where L = 35 mm, this condition required a height such that h Յ 1.75 mm. We choose h =1 mm ͑ = 0.04͒ in the present simulation, consequently the Bond number is fixed at Bd = 0.125. For this set of parameters, a typical basic flow is obtained for ⌬T = 0.1 K ͑Ma= 7.34͒.

Analytical and numerical vertical temperature profiles are shown in Fig. 4 for three different radii. The linear dependence in z of the temperature profile is clear. The first-order solution is more accurate but is not conclusive because of the strong variation in z of the zeroth-order solution. Radial velocity profiles at three different radii are presented in Fig. 5. As the Bond number is much lower than 1, the thermocapillary effects prevail. The expression of u 0 , Eq. ͑3͒, can be reduced only to the Marangoni term:

u 0 * ͑x * ,z * ͒ Ӎ -B ⌽͑x * ͒ Pr Ma 4 ͑3z *2 -2z * ͒.
This expression represents the well-known parabolic profile of a return flow with a x-dependent term ⌽͑x͒. For large values of h, comparison between experimental data due to Favre and numerical simulations are shown in Figs. 6 and7. The conditions for both simulation and experiment are h = 14.9 mm ͑Bd= 27.7͒, ⌬T = 2.9 K ͑Ma= 47255͒, moreover, for the simulation, the free surface is assumed to be insulating ͑Bi= 0͒. The analytical solution is less accurate in this case because the aspect ratio is too large ͑Ͼ0.05͒.

Figure 6 gives the experimental and numerical radial velocity as a function of the depth z. The velocity profiles are no longer parabolic because the Bond number is high, and, furthermore, gravity plays a significant role. The agreement between the numerical and experimental results indicates that the assumptions of free surface, flat and adiabatic, are relevant for the resolution of the flow.

Figure 7 presents the experimental and numerical isotherms of the basic flow in the same configuration. In the lower zone, horizontal isotherms indicates that the bottom is cold. A thermal boundary layer develops along the cold vertical wall on the right. A difference with the numerical simulation is the behavior of the isotherms near the free surface. Some experimental isotherms hatch on the free surface which proves that it is not strictly insulating contrary to the assumptions made in the simulation. It would be necessary to take into account the thermal losses by convection over the free surface, the problem being to determine a value of the heat transfer coefficient ͑or Bi͒.

B. Two-dimensional instabilities

According to the stability diagram, the primary instability corresponds to small heights of fluid, h Յ 3 mm. The basic flow, which is a return flow, destabilizes and gives birth to rolls of a circular axis surrounding the heater. For approximately ⌬T = 1 K, the first roll appears near the hot side as predicted by the theory ͓13͔ for fluids of Prandtl number PrՆ 4. Instability propagates then toward the cold edge with the appearance of new corotating rolls. The critical value of ⌬T ͑or Ma͒ for which the first roll appears is hard to determine because the rolls seems to be a part of the basic flow. Moreover, experimental thresholds reported by Favre were found using a low accuracy shadowgraphic visualization method. So values reported on the stability diagram ͑3.4 K in this case͒ are greater than the numerical ones because Favre could simply not visualize very small rolls. The number of corotating rolls depends on the Bond number and the aspect ratio of the cavity. The numerical simulation indicates that the aspect ratio influences only the number of corotating rolls present in the cavity. On the contrary, the Bond number does not change the wavelength of the instability but there is a critical value ͑Bd c Ӎ 1͒ beyond which the corotating rolls could not exist whatever the value of the Marangoni number. This can be explained by the fact that the corotating roll is a thermocapillary instability, when BdϾ 1, the thermocapillary forces being weaker than the thermogravitational force.

For the experimental conditions h =2 mm ͑Bd= 0.5͒, when ⌬T is increased above a certain secondary threshold ⌬T c , bidimensional oscillating instabilities are observed in the computations. Instabilities grow starting from a basic FIG. 6. Experimental and numerical radial velocity profiles at three different radii; h = 14.9 mm ͑Bd= 27.7͒, ᭝T = 2.9 K ͑Ma = 47 255͒.

5 mm 5 mm FIG. 7. Experimental ͑top͒ and numerical ͑bottom͒ isotherms of the basic flow for h = 14.9 mm, ᭝T = 2.9 K ͑Ma= 47 255͒ ͑20 isotherms between 0.04 and 1.54 K͒.

flow obtained by a stationary resolution of the governing equations. The disturbance appears at its early stage as sinusoidal oscillations whose amplitude grows exponentially ͑see Fig. 8͒. It can be expressed for the temperature, for example, as T͑r,,z,t͒ = T͑r,,z͒e ␤ I t cos͑␤ II t͒, ͑5͒

where ␤ I is the growth rate of instability and ␤ II represents the oscillation frequency. The growth rate value of the instability is easily obtained in a semilogarithmic graph. After this growing phase, the instability stabilizes in nonsinusoidal oscillations with constant amplitude. Figure 9 shows temporal evolution of the radial velocity at a monitoring point P ͑r = 30 mm and z =1 mm͒ over one period. The destabilization came from the first roll and consists in a periodic release of a new roll starting from the main one.

The new threshold, depending on the Bond number, could be determined precisely with a method due to Li et al. ͓24͔. Plotting ␤ I as a function of ⌬T ͑Ma͒ allows us to determine, by linear regression, the critical value of the Marangoni number ͑Ma c ͒ for which instability occurs for ␤ I = 0. In Fig. 10, we can deduce a critical temperature difference of 9.4 K ͑Ma c = 2760͒. The value of ⌬T corresponding to the appearance of hydrothermal waves in the stability diagram by Favre is 4 K ͑Ma= 1150͒, so Favre did not mention this temporal two-dimensional instability because hydrothermal waves appeared before, and experimentally these two-dimensional instabilities could not be observed. These two-dimensional oscillations were not the subject of a thorough study because Favre did not observe them. A lot of numerical publications already studied them in the past ͓14-16͔ because they do not require three-dimensional simulations. Experimentally some authors have reported them in a rectangular pool and they are often erased by three-dimensional instability.

C. Three-dimensional instabilities: Comparison with experimental data

High height of fluid: Ray-type instability

The "ray"-type instability has been simulated only in a fraction of 36°of the cavity with h = 14.9 mm. Indeed, due to the large height of fluid, we had to consider a vertical heater and periodic conditions on the lateral boundary to allow a fine mesh and reasonable computation times ͓160͑r͒ ϫ 64͑͒ ϫ 60͑z͔͒. We have chosen to model 36°of the cavity because the wave number predicted by Favre for this kind of instability was 10. Moreover, computations were performed with a cavity with a maximal angle of 45°but in this case the ray instability does not appear. For ⌬T = 5 K, isotherms on the surface obtained with a stationary computation are shown in Fig. 11. In order to obtain a similar visualization of the rays, the azimuthal temperature gradient ͑ץT / ͒ץ is computed. Moreover, an annular zone near the heater has been hatched because the shadow-graphics method was perturbed by the inner meniscus and the heater appears larger than it really is in the experimental visualization. A reconstitution to 360°of the azimuthal temperature gradient is compared to the shadow-graphic view obtained experimentally by Favre in Fig. 12. The unstationary secondary instability of the rays described by Favre corresponding to the oscillation of the structure around their initial position had not been simulated because it needed to simulate the full cavity.

Intermediate height of fluid: Petal-type instability

The petal instability has been recovered by a threedimensional simulation of a 4-mm-height layer ͑Bd= 2͒ with a horizontal heater and a cooled bottom boundary. The mesh used for the simulation is 80͑r͒ ϫ 120͑͒ ϫ 30͑z͒ and 80 000 more cells under the heater. For ⌬T =10 K ͑Ma= 11 744͒, a visualization of azimuthal temperature gradient ͑ץT / ͒ץ is compared to the shadow-graphic view obtained experimentally by Favre in Fig. 13. The experimental visualization exhibits 14 structures against 16 for the numerical one. For this value of the Marangoni number, the numerical structures are stationary contrary to the prediction of the diagram of stability in which for ⌬T Ն 6.5 K the petals rotate around the heater. The flow near the heater is quite axisymetric and bidimensional, and instability occurs at r = 20 mm, corresponding to the limit of the zone hidden by the inner meniscus. The fluid near the top surface moves faster periodically with . This produces a horizontal return flow between theses structures. In fact, these petals are similar to non-fullydeveloped rays. For this intermediate height of fluid this instability fills all the height of fluid but is unable to reach the outer wall because the bottom is cooled.

Small height of fluid

Primary instability: The co-rotative roll instability (CR).

The first observed instability is the corotative roll twodimensional instability already discussed in Sec. III B. The visualization of filled contours of the radial thermal gradient ͑ץT / ץr͒ is similar to the experimental shadow-graphic visualization of this instability as shown in Fig. 14. As predicted by Favre, the two-dimensional ͑2D͒ oscillating instability has not been observed in the 3D simulation because hydrothermal waves appeared.

Secondary instability: Hydrothermal waves (HW). The hydrothermal waves are more difficult to simulate because they need a full 3D simulation with a fine discretization along the direction. The hydrothermal waves appear in the same way as the bidimensional oscillations and stabilize after around 200 s of calculation. The same technique is used to determine the critical Marangoni number value for which they appear. The initialization field of the computation is deduced from a stationary result with the same conditions which is normally an axisymetric flow.

In the experimental configuration using a 120͑r͒ ϫ 240͑͒ ϫ 15͑z͒ plus 48 990 meshes under the heater, hydrothermal waves were found for h =1 mm ͑Bd= 0.127͒ and ⌬T =5 K ͑Ma= 367͒. Figure 15 shows three-dimensional isosurfaces of azimuthal thermal gradient ͑ץT / .͒ץ These surface waves look like curved branches formed in a wave source and disappearing in a wave sink. These waves propagate perpendicularly to the isotherm as in Favre et al. ͓2͔. Two wave sources and two wave sinks were found. Hydrothermal waves obtained for other ⌬T and different height were always with multiple sources and sinks.

Three simulations in this configuration for ⌬T = 3, 5 and 7 K allowed us to estimate the growth rate of this 3D hydrothermal wave as a function of ⌬T. The critical temperature difference for wave apparition is ⌬T c = 3.8 K ͑Ma= 279͒. From the diagram of stability of Favre et al. ͓2͔ for a liquid height equal to 1 mm, the hydrothermal waves appear at the critical difference of temperature of 3 K ͑Ma= 220͒ in agree-ment with the simulations. Note that Shi and Imaishi ͓27͔ found a similar critical threshold for the apparition of hydrothermal waves ⌬T c = 5.9 K for a configuration with Pr= 6.7, h = 1 mm, R 1 = 20 mm, R 2 = 40 mm with an adiabatic bottom and where the fluid is heated by the outer wall and cooled by the inner wall.

The spatiotemporal experimental and numerical diagrams of these hydrothermal waves are shown in Fig. 16. The numerical diagram was built with the fluctuating temperature TЈ, defined by Eq. ͑2͒, and computed on a circle at r = 20 mm on the free surface. The experimental diagram exhibits one wave source and one wave sink whereas the numerical diagram shows two sources and two sinks for the hydrothermal waves. The hydrothermal waves are a 3D extension of the corotating rolls as an azimuthally thermal "strength" where the ⌬T increases.

Simulations of a higher layer of fluid corresponding to the experimental setup ͑h =2 mm͒ do not give proper results. In fact, the initial stationary computation does not give axisymetric results. Therefore the study of appearance of HWs were impossible. We found that this problem does not occur when the fluid is not cooled by the lower wall; in this configuration the radial extension of the instability is higher.

The main difference between simulation and experiment is the adiabatic free surface ͑Bi= 0͒ which may change the number of branches but not the temperature threshold for the apparition of the waves. A fine numerical simulation of thermoconvective instabilities of a fluid with an intermediate Prandtl number requires important resources. The buoyant-thermocapillary instabilities in a cavity filled with a fluid of Prandtl number of 10.3 heated on the surface in the center and cooled at the bottom and lateral boundary were simulated and compared with experimental results with success. The two-dimensional and three-dimensional simulations of the basic flow correspond to the experimental results of Favre ͓2͔ and of Castem 2000 simulations ͓3͔. For the small aspect ratio of the cavity, the analytical solution of Garnier ͓9͔ was extended to the specific case of a cooled lower wall and was verified by the simulation. The different stationary instabilities observed by Favre have been studied and successfully simulated. The ray instability has been described for a high height of fluid. For an intermediate height of fluid petal instability has been simulated. We have shown that for small height of fluid the threshold of two-dimensional corotating instabilities is greater than that of the three-dimensional hydrothermal waves. It explains why experimentally no unstationary twodimensional instability was reported in this case by Favre ͓2͔. Like in experimental observations, hydrothermal waves appear only when the dynamic Bond number is less than 1, i.e., when thermocapillary effects are predominant. Hydrothermal waves were more difficult to simulate, and we had to change the configuration of the heater to obtain valuable results. It can be explained by the fact that the inner and outer meniscus were not simulated. Moreover, the experimental heating boundary conditions are not strictly an imposed temperature.

  FIG. 1. Diagram of stability due to Favre ͓2͔ ͑see text for details͒.

  1 W m -1 K -1

  FIG. 2. Shadow-graphic views of stationary instabilities: ͑a͒ corotating rolls ͑CRs͒, ͑b͒ stationary petals ͑SPs͒, and ͑c͒ stationary rays ͑SRs͒ obtained by Favre ͓2͔.

FIG. 4 .FIG

 4 FIG. 4. Nondimensional axial temperature profiles at three different radii ͑25, 30, and 38 mm͒ and Bi= 5; h =1 mm ͑Bd= 0.125͒, ᭝T = 0.1 K ͑Ma= 7.34͒.

  FIG. 8. Typical temporal evolution of the corotating instability.

  FIG. 10. Growing rate of the bidimensional instability ͑␤ I ͒ for different applied thermal difference ⌬T.

FIG. 12 .

 12 FIG.12. Numerical ͑⌬T =5 K, Ma=81 475͒ and experimental ͑⌬T = 10 K, Ma= 162 950͒ visualization of the stationary ray-type instability for h = 14.9 mm ͑Bd= 27.7͒, the hatched zone corresponds to the zone hidden by the inner meniscus in the experimental visualization.

FIG. 14 .

 14 FIG. 14. Stationary corotative roll, experimental ⌬T = 3.2 K, ͑Ma= 761͒ and h = 1.8 mm ͑Bd= 0.41͒ shadow-graphic and numerical ⌬T =7 K ͑Ma= 2055͒ and h =2 mm ͑Bd= 0.5͒ visualization, scale between 0 and 100 K m -1 .

  FIG. 16. Spatiotemporal diagrams: ͑a͒ experimental for h =2 mm ͑Bd= 0.51͒ and an unknown ⌬T, and ͑b͒ numerical for h =1 mm ͑Bd= 0.127͒ for ⌬T =5 K ͑Ma= 367͒.

TABLE I .

 I Physical properties of the fluid.

TABLE II .

 II Grid convergence for the two-dimensional case.

	Number
	of meshes	u͑25; 14, 9͒	u͑30; 14, 9͒	u͑38; 14, 9͒		% mesh	% mesh
	͑r ϫ z͒	mm s -1	mm s -1	mm s -1	Pe Max	PeՅ 1	PeՆ 10	Flux ͑W͒
	7 636	5.2	2.6	1	35.8			0.2197
	50 700	6.6	4	1.1	17.3	75	1.5	0.2237
	110 000	7.3	4.8	1.1	15	77	0.3	0.2258
	126 000	7.6	5	1.9	13	78	0.2	0.2265

TABLE III .

 III Grid convergence for the 3D case with two meshes: A is 80͑r͒, 120͑͒, 20͑z͒, and B is 120͑r͒, 180͑͒, 18͑z͒. Around 50 000 meshes are used in both case A and B to discretize the central zone under the heater.

	tal case of Favre which considers a cooled bottom boundary.
			u͑25; 2͒	u͑30; 2͒	u͑38; 2͒			% mesh	% mesh
	T ͑K͒ Mesh	mm s -1	mm s -1	mm s -1	Pe Max Pe average	PeՅ 1	PeՆ 10	Flux ͑W͒
	1	2D	1.60	1.25	0.7					0.0674
		3D A	1.55	1.2	0.68	23	1.9	48	2.9	0.0662
		3D B	1.55	1.15	0.65	21	1.39	48	1.17	0.0651
	2.9	2D	3	2.5	1.75					0.227
		3D A	2.9	2.4	1.75	48	4.64	24	17	0.222
		3D B	2.8	2.4	1.75	43	3.3	20.6	8.9	0.227
	5	2D	2.5	4	2.5					0.420
		3D A	3.5	3.25	2.5	68	6.73	20	24	0.413
		3D B	3.4	2.75	2.2	98	6.4	8.9	22	0.457

TABLE IV .

 IV Adimensionnal scale. Temperature ͑T -T c ͒ / ͑T H -T c ͒

	COMPARISON BETWEEN NUMERICAL AND EXPERIMENTAL …	PHYSICAL REVIEW E 76, 066307 ͑2007͒
	Radial length ͑x͒	L
	Height ͑z͒	h
	Velocity ͑u , w͒	/ h
	Pressure	2 0 L / h 3
	Density	0
		066307-5


PhD student Julien Ogor
email: julien.ogor@ensta-bretagne.org

Benoit Zerr
email: benoit.zerr@ensta-bretagne.fr

Towards the automation of sand dune detection in the bathymetry

Introduction

The advances in echo sounder (MultiBeam EchoSounder) and positioning (GPS) technologies of the past two decades enabled us to survey larger areas, faster and more accurately. The collected data show how diverse underwater sand dunes are in terms of shape, spatial organization, dimensions, dynamics, etc. To deal with this large amount of bathymetric data, manual processing appears outdated. By consequent, tools allowing to quantitatively describe their morphology and dynamics have been developed [1,2,3]. It mainly enables to objectively characterize the sand dunes and banks. Unlike the existing tools, we aim at describing each dune independently with a set of morphometric parameters. With this goal of dune characterization, the first step consists in detecting the dunes in the bathymetry. As a preliminary approach, we propose a detection technique inspired from geomorphometry.

Geomorphometric approach

Results

Perspectives

[1] Dorst, L.L. 2004. Geodetic deformation analysis: new method for the estimation of seabed dynamics; Proc. MARID Workshop, Enschede, 1-2 April 2004.

Datasets References

Many thanks are due to the FPS Economy Continental Shelf Service and SHOM for giving the author a full access to their data. Also, the author is grateful to the DGA and Region Bretagne for their financial support..

[ 2 ]

 2 Lisimenka A. & Rudowski S. 2013. Bedform characterization in river channel through 2D spectral analysis. Proc. MARID Workshop, Bruges, 15-16 April 2013. [3] Knaapen M.A.F, van B. Henegouw C. N. & Hu Y. Y. 2005. Quantifying bedform migration using multi-beam sonar. Geo-Mar. Lett., vol. 25, no. 5: 306-314, Nov. 2005. [4] Grosse, P. & al. 2012. Systematic morphometric characterization of volcanic edifices using digital elevation models. In Geomorphology, 136 (1): 114-131. [5] Graff, L. H. 1992. Automated classification of basic-level terrain features in digital elevation models. DTIC Document. [6] d'Oleire-Oltmanns S., EisankC., Dragut L. & Blaschke T. 2013.An Object-Based Workflow to Extract Landforms at Multiple Scales From Two Distinct Data Types. In IEEE Geoscience and Remote Sensing Letters, 10(4): 947-95

Acknowledgements

Marine & River Dune Dynamics Conference V, 4 th -6 th , April 2016 . North Wales, UK.

Dune Delineation

Among the large number of local characteristics commonly used in geomorphometry to detect features at the Earth's surface (e.g. volcanoes [4], mountains [5], drumlins [6]), the slope and minimum curvature seem the most appropriate to our case : dunes.

Local Characteristics Selection Dune Delineation

To conduct this research, SHOM (French Navy Hydrographic Office) and the FPS (Federal Public Service) Economy gave us access to their dune databases. These two datasets are complementary as they were not collected in the same prospect.

SHOM Dataset

Long lines run to explore the continental shelf (Brittany, straits of Dover) or boxes surveyed over interesting structures Depth ranging from 30 to 200 m Large to giant dunes ( below 1m high to over 20 m high) Data acquired during annual campaign

FPS Economy Dataset

11 small areas spread around the Belgian continental shelf surveyed to monitor the impact of aggregate extraction Depth ranging from 5 to 40 m Few meter-high dunes with sand waves superimposed. Surveys completed several times every year since 1999.

These datasets are perfect to appreciate the wide variety of dunes in terms of morphology and dynamics.

The region growing algorithm is divided into two phases: 1. Define the seed regions (Crest Regions)

2. Find a criterion to stop the region growth

The preliminary results are promising and confirm that the choice of designing a geomorphometric approach to automatically detect dunes is realistic. Indeed, the algorithm was tried out on a dataset for which we also have manually-drawn crest lines. 248 dunes are detected by the algorithm against 252 with the classical, manual approach.

The minimum curvature is very effective to highlight the crest locations. To base the growth stage on the slope is too simple as such. Another metrics, criterion inspired from the slope is to be found.

It is essential that the crest regions are perfectly extracted prior to the region growth. Furthermore, the region growth is very challenging as the dune limits are sometimes not obvious to a human expert.

In addition, the dune diversity of shapes makes impossible the use of crisp thresholds appropriate to all the dunes.

Conclusion

Things to be considered in order to design a fullyautomated dune detection algorithm:

Bathymetry

Crest Regions Dunes

Depth (m)

2013.

Observations:

Most dune crests are correctly detected Some crests are partly detected or/and split into few regions (Frame 1)

A crest region can gather several dune crests (Frame 2)

Observations:

Slope varies from a dune to another or even from a flank to the other: