## Abstract

The existence of an oscillatory instability in the Bénard–Marangoni phenomenon for a viscoelastic Maxwell's fluid is explored. We consider a fluid that is bounded above by a free deformable surface and below by an impermeable bottom. The fluid is subject to a temperature gradient, inducing instabilities. We show that due to balance between viscous dissipation and energy injection from thermal gradients, a long-wave oscillatory instability develops. In the weak nonlinear regime, it is governed by the Korteweg–de Vries equation. Stable nonlinear structures such as solitons are thus predicted. The specific influence of viscoelasticity on the dynamics is discussed and shown to affect the amplitude of the soliton, pointing out the possible existence of depression waves in this case. Experimental feasibility is examined leading to the conclusion that for realistic fluids, depression waves should be more easily seen in the Bénard–Marangoni system.

## 1. Introduction

The study of solitons in physical systems is intimately connected to the Korteweg–de Vries equation (KdV), whose integrability has been established in Gardner *et al*. (1967). KdV is usually derived as a result of asymptotic expansions in conservative systems, as in the case of the long waves in shallow water approximation to the free surface Euler's equations in hydrodynamics. In this case, solitons appear most commonly as elevation waves, depression waves being rather rare. In this work, we will point out the possibility of creating depression solitons on the surface of a fluid. However, the system we will study is not a genuinely conservative system. We will consider a viscous fluid that is heated from below, allowing a balance at a certain critical number so as to ensure the existence of solitons. Moreover, we will consider viscoelastic effects, which will be shown to be connected to amplitude inversion. In this context, wave propagation on the surface of a fluid is connected to the existence of purely oscillatory instabilities of the system, in the linear limit, with vanishing wavenumber.

We will address in this work the question of the evolution of a free surface of a viscoelastic Maxwell's fluid heated from below, under stress-free mechanical conditions and certain thermal boundary conditions. The influence of the temperature field comes in through two different mechanisms: buoyancy, meaning the temperature dependence of the density, in which case the system is called Rayleigh–Bénard, or the Marangoni effect, consisting of the temperature dependence of the surface tension, which goes under the name Bénard–Marangoni system. Here, we will consider both effects jointly. Moreover, we will be interested in the surface perturbations of a fluid in a purely conductive state, that is, in the absence of convection. Therefore, the non-dimensional numbers connected to the onset of fluid motion, the Rayleigh number *Ra* and the Marangoni number *Mg*, will be smaller than their critical values for convection.

In the case of a Newtonian fluid, the existence of oscillatory instabilities of the free surface is a well-known result of the theory of hydrodynamic stability. We will concentrate here on the case in which the temperature is fixed at the bottom and the temperature gradient is fixed at the surface. At the bottom, no-stress boundary conditions are used. In this case, for the joint Rayleigh–Bénard–Marangoni system, the linear analysis shows (Benguria & Depassier 1989; Chu & Velarde 1989) that for *Ra*−(5/2)*Mg*=30 a long-wave oscillatory instability of the motionless state sets in. This instability is not connected to convection. Rather, at the onset point, what we have is a balance between viscous dissipation and energy injection from the thermal gradient, sustaining waves on the surface. In the weak nonlinear limit, this instability is characterized by existence of solitons propagating at the free surface (Alfaro & Depassier 1989; Kurcbart *et al*. 1990; Chu & Velarde 1991; Kraenkel *et al*. 1992*a*). Near criticality and taking higher order terms into account, the theory predicts the existence of stable dissipative–dispersive solitary waves (Marchant & Smith 1998; Depassier 2006), sometimes also called dissipative solitons (Akhmediev & Ankiewicz 2005). In fact, it is the temperature gradient that balances viscous dissipation and sustains the waves. The actual existence of these objects has been established experimentally for the case of the surface tension-dominated case in Weidman *et al*. (1992) and Linde *et al*. (2001).

In Ramkissoon *et al*. (2006), the linear stability analysis for the analogous situation with a viscoelastic Jeffrey fluid was performed. Results for the Maxwell fluid limit were deduced as a special case. It was shown that the critical value for the onset of the oscillatory instability is not affected by viscoelasticity. Therefore, in such fluids, the same long-wave instability also exists. In this work, we extend the previous linear result on the Maxwell fluid to the weakly nonlinear regime in order to unveil the viscoelastic effects on the surface waves that result from the instability. It will turn out that viscoelasticity introduces new dispersion into the system, which will alter the evolution of surface waves. For a given amplitude, the corresponding solitary wave will be seen to be longer for a Maxwellian fluid than for a Newtonian one. Considering one-soliton solutions, for a strong enough viscoelastic effect, a wave of depression, instead of elevation, will be possible. We note that depression surface waves are a rare phenomenon, having been observed in very thin layers of mercury (Falcon *et al*. 2002). Our results imply that exploration of viscoelastic effects should lead to similar situations in less thin fluids.

We should point out that previous work has appeared in connection with viscoelastic fluids acted upon by thermal gradients. For instance, Lebon *et al*. (1994), Parmentier *et al*. (2000) and Li & Khayat (2005) considered convective instabilities that occur in such systems, in particular due to buoyancy and surface tension effects. We note, however, that the instability addressed in the present paper is of a different kind, connected to deformable free surface effects, whose nonlinear stage has not been considered earlier, and setting in for thermal gradients much below the ones necessary to produce convection.

This paper is organized as follows: in §2, we give the physical and mathematical model explicitly and introduce non-dimensional variables, rendering the equations amenable to a perturbative treatment; in §3, a weakly nonlinear long-wave expansion is performed and the equations for the free surface elevation are found; in §4, the theoretical results are analysed; in §5, we discuss the experimental feasibility; and in §6, we make our final remarks.

## 2. The Rayleigh–Bénard–Marangoni system: mathematical formulation

We consider a layer of viscoelastic Maxwell fluid initially at rest in a Cartesian frame (*x*, *y*, *z*), so that (*Oz*) is the upward vertical direction of unitary vector . We assume no dependence on *y*, allowing us to consider a two-dimensional system in the *x*–*z* plane. The fluid domain is contained between a bottom at *z*=0 (flat and perfect thermal conductor) and a upper free surface at in contact with a passive gas, which exerts upon it a constant pressure *p*_{a}. The undisturbed state *z*=*d* is acted on by a gravitational field . The fluid is described by the Navier–Stokes equation in the Oberbeck–Boussinesq approximation (Kraenkel *et al*. 1992*b*; Kozhoukharova & Rozé 1999)(2.1)(2.2)(2.3)(2.4)(2.5)and the constitutive equation for the upper convected Maxwell model(2.6)In equations (2.1)–(2.6), is the convective derivative; is the velocity; *p* is the pressure; *T* is the temperature; *ρ* is the density; *ζ* is the surface tension; and is the extra-stress tensor (Astarita & Marruci 1974; Bird *et al*. 1987; Joseph 1990). In the constitutive equation, D/D*t* is the upper convected derivative defined by(2.7)where superscript *t* means transposition. The fluid properties, i.e. its viscosity *μ*, thermal diffusivity *κ*, the coefficient of thermal expansion *α*, the rate of change of the surface tension with temperature *γ* and the relaxation time *λ*, are constant. The quantities *ρ*_{0}, *T*_{0} and *ζ*_{0} are reference values.

In order to obtain the dimensionless form of the equations, we introduce dimensionless quantities indicated by *x*′, *z*′, … and defined byFurthermore, we introduce the following dimensionless numbers:*σ* is the Prandtl number; *Ga* is the Galileo number; *Ca* is the capillarity number; *Ra* is the Rayleigh number; *Mg* is the Marangoni number; and *L* is the dimensionless relaxation time. *Ra* and *Mg* are two control parameters and *L* can be viewed as a Deborah number.

Substituting these into (2.1)–(2.6) and dropping the primes, we obtain(2.8)(2.9)(2.10)(2.11)(2.12)(2.13)(2.14)and the equations of state(2.15)

The static solutions of equations (2.8)–(2.15) are(2.16)

(2.17)

In order to write the boundary conditions, we need to introduce the unit normal vector ** n** and a tangential vector

**. They are defined by(2.18)where . We now complete the fundamental equations (2.8)–(2.14) with appropriate kinematic, dynamical and thermal boundary conditions (Wehausen & Laitone 1960; Batchelor 1967). The classical kinematic boundary condition for the free surface**

*t**S*(

*x*,

*z*,

*t*)parametrized as yields toWe shall assume that the fluid is bounded below by a perfect thermal conductor and subject to a fixed thermal gradient above. We consider

*T*constant on

*z*=0 and we fix the heat flux through the fluid, so that we have on where

*F*is the constant normal heat flux and

*k*is the thermal conductivity. Now, non-dimensionalizing and adopting Δ

*T*=

*Fd*/

*k*as unit of temperature, we obtainThe stress condition at the free surface (Laplace law) isSo, explicitly written, the boundary conditions on read(2.19)(2.20)(2.21)(2.22)At the bottom

*z*=0, we assume a stress-free boundary condition (otherwise oscillatory instabilities do not exist; Benguria & Depassier 1989; Kurcbart

*et al*. 1990; Ramkissoon

*et al*. 2006)(2.23)

## 3. Multiple scales expansion approach: long-wave equation

Let us now look for the nonlinear evolution of long-wave and small amplitude perturbations to the static solutions (2.16) and (2.17), so we consider(3.1)(3.2)In order to determine *u*″, *w*″, *p*″, and *θ*, we will use the very well-known multiple scales or reductive perturbation method (Taniuti & Wei 1968; Jeffrey & Kawahara 1982). We introduce appropriate slow space and time variables needed to describe long waves in *x* asymptotically in the time *t*. They are(3.3)where the small parameter *ϵ* characterizes the amplitude of the surface perturbation and *c*, its velocity, must be determined as a compatibility condition. In this work, we will consider *ζ*_{0}=0 that implies 1/*Ca*=0. Introducing now the expansions (easily justified using (3.3))(3.4)(3.5)(3.6)(3.7)(3.8)(3.9)(3.10)(3.11)where all quantities are dimensionless functions of the variables *ξ*, *z* and *τ*, we can obtain an order by order solution to the problem.

At order *ϵ*^{0}, the solution is given by(3.12)At order *ϵ*^{1}, we find(3.13)(3.14)with *f*(*ξ*, *τ*) an arbitrary function. At order *ϵ*^{2}, we find(3.15)(3.16)(3.17)(3.18)(3.19)(3.20)(3.21)(3.22)with *g*(*ξ*, *τ*) an arbitrary function. To obtain equations (3.17)–(3.22), we have used the relation(3.23)which comes from the solvability condition at this order. At the next order, the solution is given by(3.24)(3.25)(3.26)(3.27)(3.28)(3.29)(3.30)(3.31)with *m*(*ξ*, *τ*) an arbitrary function. The solubility condition at this order is(3.32)Finally, at order *ϵ*^{4}, we find(3.33)(3.34)and an equation for the derivative in *z* of *τ*_{3,xz}(3.35)

Now integrating (3.35) in *z* (using the expressions for *u*_{2}, *p*_{4} and *τ*_{1,xx}), we find(3.36)

The requirement of compatibility of equations (3.36) and (3.31) provides as evolution equation for the function *f*(*ξ*, *τ*) the KdV equation(3.37)

## 4. Discussion of the theoretical results

As *η*_{0}=*f*/*c*, equation (3.37) is also the equation for the surface displacement. We should recall that *Ra* and *Mg* are not independent, being related by equation (3.32). The above evolution equation includes nonlinearity and dispersion, represented, respectively, by the second and third terms in the equation. It can be transformed to the standard KdV equation(4.1)by a simple rescaling of the variables, namely(4.2)where(4.3)

(4.4)

The soliton solution to this equation reads(4.5)where *q* is arbitrary.

*L* is positive and contributes to diminishing the value of *Q*, which is the adimensional dispersive length scale, connected to the wavenumber *qQ*. Therefore, the smaller the *Q* the longer the wave.

Interestingly enough, a large positive value of *L* may change the sign of *Q*. If this happens, instead of an elevation wave one should see a depression wave travelling in the opposite direction. If *Mg*=0 and thus *Ra*=30, the critical value of *L* for the elevation–depression inversion is found to be(4.6)A further interesting case comes in when we consider Marangoni dominated flows with small Galileo number. In this case, the effects of gravity may be neglected. We thus can have *Ra*→0, *Ga*→0, while *Mg* (the Marangoni number) remains finite. This is typically true for sufficiently thin fluid layers, a situation that can best be tested in experiments, as *Ga*∼*d*^{3}, *Ra*∼*d*^{3} but *Mg*∼*d*. With these assumptions, we have *Mg*=−12, *c*^{2}=12*σ*, *M*=1 and . If we define a critical *L* by *Q*=0, it obeys(4.7)If *L*>*L*_{c} then we will have a depression wave instead of the usual elevation wave, which would be the only allowed case if no viscoelastic effect were considered. It follows that for Maxwellian fluids (where the viscoelastic effects are due to relaxation) that the inversion of elevation to depression wave can in principle be observed.

Integrable evolution is, in this case, a first-order property. This is in contrast to the case of Euler equations, and for conservative systems in general. It is known that (Fokas & Liu 1996; Dullin *et al*. 2003) the next order correction to the KdV equation in conservative systems has the form(4.8)which is asymptotically integrable (Kodama 1987), in the sense that there is a near-identity transformation taking it to an integrable equation, correct up to order (*ϵ*^{2}) and reducing the correction to the next commuting flow in the KdV hierarchy. Obstacles to integrability appear only at the next order. In our case, the situation is different. Although the lowest order of the long-wave expansion results in the usual KdV equation, the corrections appear at (*ϵ*) and have the general form (Aspe & Depassier 1990):(4.9)

This equation is not related to an integrable one by a near-identity transformation and represents the breaking of asymptotic integrability at a lower order than usual. This comes from the dissipative nature of the system, where the balance between viscous dissipation and energy injection through heating is only exact at the lowest order. The effects of the corrections will thus be typically of a non-integrable kind. In this case, this results in amplitude selection, an effect characteristic of dissipative solitons (Akhmediev & Ankiewicz 2005).

## 5. Experimental feasibility

In this section, we discuss on the experimental possibilities to observe surface *depression* solitons. Let us introduce the ratio *Γ* between the Marangoni number *Mg* and the Rayleigh number *Ra*Values of *Γ*>1 mean that the Bénard–Marangoni phenomenon is the dominant one, and values of *Γ*<1 mean that the Rayleigh–Bénard phenomenon dominates. With the definitions(5.1)we have(5.2)The value *Γ*=1 determines the critical depth *d*_{c}. Two very common Maxwell's fluids are *wormlike micelles* and *telechelic polymers* solutions. For these solutions, values of *ρ*_{0}, *γ*, *α* (and also *κ*) are very near of those for water. So, using figures from table 1, we obtain(5.3)Now, with this value for *d*_{c} and the integrability requirements *Ra*=30 or *Mg*=−12, we choose appropriate values of *μ* for Maxwell's solutions and realistic Δ*T* to be in the Rayleigh–Bénard or in the Bénard–Marangoni conditions. This is summarized in table 1. In the upper row, we have listed the values of physical parameters for Rayleigh–Bénard and in the lower row those for Bénard–Marangoni.

The possibility or not to observe a *depression* soliton depends on the values of *λ* determined by the expressions (4.6) and (4.7) through(5.4)Both phenomena bring to the same (approx.) *L*_{c}∼0.4, so we have

in the Rayleigh–Bénard case:

*λ*>10^{4}sin the Bénard–Marangoni case:

*λ*>10^{2}s.

Hence the Rayleigh–Bénard case necessitates an enormous relaxation time (*λ*∼2.7 years), which is non-realistic. Thus, the Bénard–Marangoni case is the appropriate one to observe depression surface solitons (*λ*∼1.6 min).

## 6. Final comments

In summary, we have discussed the existence of solitary wave excitations on the surface of a heated viscoelastic fluid (Maxwell fluid). We obtained conditions for the existence of such waves and discussed the relevance of the viscoelastic effects, pointing out the possibility of having depression waves instead of elevation ones in such systems. The system under inspection behaves as integrable only at the lowest order, and asymptotic integrability breaks at the next order, in contrast to the usual conservative case. A feasibility analysis shows that for acceptable values of the physical parameters the existence of the solitary waves is more likely for the case of the Bénard–Marangoni system.

## Acknowledgments

M.A.M and D.C. wish to thank the Instituto de Física Téorica-UNESP for hospitality and partial support. R.A.K. thanks CNPq (Brazil) for partial support. The authors also thank S. Mora from Laboratoire des Colloïdes, Verres et Nanomatériaux-UM2 (France) for discussions about experimental feasibility.

## Footnotes

- Received May 27, 2008.
- Accepted August 13, 2008.

- © 2008 The Royal Society