STABILIZING INTERPLAY BETWEEN THERMODIFFUSION AND VISCOELASTICITY IN A CLOSED-LOOP THERMOSYPHON

Viscoelastic fluids represent a major challenge both from an engineering and from a mathematical point of view. Recently, we have shown that viscoelasticity induces chaos in closed-loop thermosyphons. This induced behavior might interfere with the engineering choice of using a specific fluid. In this work we show that the addition of a solute to the fluid can, under some conditions, stabilize the system due to thermodiffusion (also known as the Soret effect). Unexpectedly, the role of viscoelasticity is opposite to the case of single-element fluids, where it (generically) induces chaos. Our results are derived by combining analytical results based on the projection of the dynamics on an inertial manifold as well as numerical simulations characterized by the calculation of Lyapunov exponents.

1. Introduction.Thermosyphons are special heat exchange devices that function autonomously by exploiting the effect of natural convection (due to temperature differences) on the motion of a fluid.
Many oils and coolants used in thermosyphons are composed of long-chain polymers, having viscoelastic properties.Unlike Newtonian fluids, viscoelastic fluids share some properties with elastic solids.In particular, they have both viscous and elastic characteristics.Viscous materials, like water, cannot resist shear and deform linearly when a shear stress is applied [1].Elastic materials, like solids, strain instantaneously when stretched and, unless the plastic limit is reached, they cannot flow.Elasticity is the result of bond stretching along crystallographic planes in an ordered solid, whereas viscosity is the result of the diffusion of atoms or molecules inside an amorphous material.Depending on the change of strain rate versus stress inside a material, the viscosity can be categorized as having a linear or nonlinear response.When a fluid exhibits a linear relation between stress and instantaneous strain rate, it is called, generically, a Newtonian fluid.In this case the stress is linearly proportional to the strain rate.If the material exhibits a nonlinear response to the strain rate, it is called a non-Newtonian fluid [6].Besides the linear response to stress, viscoelastic materials possess elastic restoring forces like solids do.In summary, viscoelastic materials have elements of both, viscous as well as elastic properties and, as such, exhibit time dependent strain.This combined effect of both liquids and solids can be quantified experimentally by means of constant stress experiments (creep) or relaxation after stress removal [6].
In a recent work [32], we showed that viscoelasticity can induce chaotic behaviors related to the memory effects intrinsic to this sort of material.Thus, the stronger the influence of viscoelasticity, the faster is the route to chaos.This situation poses an appealing question from an engineering perspective: under which conditions or mechanisms can we remove this chaotic oscillation?This question arises naturally as, in practice, thermosyphons are expected to operate under stable conditions (namely, the thermosyphon will be more efficient if the velocity sustained inside it is constant).
One possibility aimed to control this chaotic behavior is to explore their response to non-equilibrium conditions like thermodiffusion (also known as thermophoresis or the Soret effect), besides the intrinsic properties of the material.Thermodiffusion refers to the effect caused by a temperature gradient [4] in a mixture of two or more diffusing substances.The term "Soret effect" normally means thermodiffusion in liquids [12].
The Soret effect is a very important component of the study of any physical experiment that incorporates thermodiffusion.It gives rise to interaction between the thermal and solute gradients even when the fluid is at rest [19].Inside a thermosyphon, because of the temperature gradients, the Soret effect induces solute concentration gradients, thus triggering natural convection inside the loop.
The diffusion induced by the temperature gradients and solute gradients of the viscoelastic fluid by the Soret effect is the main focus of our study.The study of the spatiotemporal phenomena in viscoelastic fluids emerges from the integrated observations of velocity, temperature and solute concentration [12].Interestingly, the effect involves a new field for the solute concentration (S, see Eq. (10)) that is doubly coupled with temperature through the equations for the velocity and solute concentration respectively.
Our contributions in this paper are: • To introduce a system of equations (10) governing a closed loop thermosyphon model of a viscoelastic fluid with the Soret effect which, although is a generalization of the previous models [32,9,16,15,14,17,22,11,28,24], increases the order of the time derivatives in velocity giving rise to different routes to chaos.• To generalize previous results for the dynamics on the inertial manifold in the case of having an additional partial differential equation for the solute concentration, coupled with both the temperature and the velocity inside the loop.• To provide a detailed numerical analysis of the behavior of the velocity, temperature and solute concentration which includes a thorough study of the various behaviors of the system for different values of viscoelastic fluid and the Soret coefficient.• To assess the interplay between the viscoelastic memory effect and the Soret effect on the system.This paper contains four sections.The first section provides a brief introduction to the problem, explaining briefly the underlying Physics of a thermosyphon with viscoelastic fluids and including the Soret effect.The second section illustrates the mathematical formulation of the model of the system, deriving a nonlinear system of equations and the physical meaning of the relevant parameters.The third section contains the numerical approximation of the nonlinear system of equations.The fourth section deals with the numerical analysis of the reduced model and discussion of the results using Lyapunov exponents.Finally, we conclude by discussing the main effects of viscoelasticity under non-equilibrium conditions and summarising our findings along with some insights for future research.
2. Mathematical formulation of the model.We consider a thermosyphon model in which the confined fluid is viscoelastic (for instance, consider the behavior of jelly or toothpaste).This has some a-priori interesting peculiarities that could affect the dynamics with respect to a Newtonian fluid (like water).On the one hand, the dynamics has memory (see below) as its behavior depends on the whole past history.On the other hand, for small perturbations, the fluid behaves like an elastic solid with a characteristic frequency of resonance that, eventually, could be relevant.
The simplest approach to viscoelasticity comes from the so-called Maxwell constitutive equation [21,3].Although this model is a great simplification, it has been proven valid even for complex fluids as blood in which red cells change its behavior depending on their concentration or even the geometry of the vessel [27].In this model, both Newton's law of viscosity and Hooke's law of elasticity are generalized and complemented through an evolution equation for the stress tensor, Σ.The stress tensor comes into play in the equation for the conservation of momentum: where p is the hydrostatic pressure, ρ the density of the material and g is the acceleration due to gravity (g = 9.81 m/s 2 ).The latter equation is supplemented with the hypothesis of incompressibility (accurate enough for liquids), referred to as the continuity equation: For a Maxwellian fluid in a thin section thermosyphon, the stress tensor, Σ, is reduced to just one independent component, say σ, that evolves according to: where µ is the fluid viscosity, E the Young's modulus and γ is the only non-zero component of the shear strain rate tensor, Γ (or rate at which the fluid deforms) defined as Under stationary flow, equation (3) reduces to Newton's law and, consequently, Eq. ( 1) reduces to the celebrated Navier-Stokes equation.For short times where impulsive behavior from rest can be expected and equation (3) reduces to Hooke's law of elasticity.
In order to describe de physical mechanisms introduced in Section 1, we will formulate a mathematical model for the following observables: the velocity v, the distribution of the temperature T of the fluid and the solute concentration S into the loop, the latter governed by transport equations.Specifically, where J T,S represent the heat and solute fluxes (or currents), and K T,S , external temperature or solute sources, respectively.Typically, those terms are specified using constitutive equations that are linear combinations of the gradients of other variables.For instance, Fourier's law of heat conduction establishes that J T = −ν∇T , so Eq. ( 4) provides the well-known Laplacian ν∇ 2 T , etc.
In order to proceed, we follow the same procedure as [29] (see Fig. 1 for a rationale of the method and some basic notation): (i) Project equations ( 1), ( 2), ( 4) and ( 5) along the tangent of the thermosyphon pipe (whose coordinate we call x hereafter).(ii) Assuming that the section, A, of the thermosyphon is small enough, we can assume that differences in the cross direction of the pipe are negligible so averaging does not produce a significant loss of accuracy.(iii) As the averaged velocity is the same throughout the pipe (because the fluid is incompressible).
Figure 1.A section of the thermosyphon: we project Eqs. ( 1), ( 2), ( 4) and ( 5) along the tangent, ŝ and average them over a section of diameter A 1. As shown, we denote x the coordinate along the thermosyphon pipe.
Before following these steps, we eliminate the stress component σ in the following way: we first differentiate Eq. ( 1) with respect to time and multiply it by µ/E.Then add the former Eq. ( 1) and use Eq. ( 3) to find Following Refs.[29,18], we assume that the non-linear terms and the viscous term, ∇ 2 v, can be expressed in terms of the so-called wall law force, F, that is a function of the average velocity (see below for details).
Because of the incompressibility hypothesis, Eq. ( 2), the averaged velocity v is the same throughout the thermosyphon, thus, it only depends on time, t, but not on the spatial coordinate along the pipe, x.Also due to this projection and averaging steps, temperature, T , and solute concentration, S, will only be functions of time, t, and position along the loop, x.
Projecting along the direction of the pipe, we arrive at: where F s is the projection of the wall law force along the pipe (see below for details about this function).Velocity in Eqs. ( 1)-( 3) is coupled with temperature and diffusion through the body force ρg, where ρ depends on density (hotter fluids are less dense) and solute concentration (higher concentrations of solute imply higher densities).Mathematically, one can assume that there is a linear relationship between density and temperature.Likewise, in this work we assume also a linear relationship between density and solute concentration.So, in the end where α T and α S are dilatation and thermophoretic coefficients.Here, (as in Refs.[29,18]) we apply the Boussinesq approximation, that states that the density, when it is proportional to the dynamical terms d 2 v/dt 2 and dv/dt, is approximately constant and equal to ρ 0 , and also that inertial contributions from temperature and solute concentration are negligible.Thus, integrating Eq. ( 7) along the length of the pipe L (so the integrated pressure gradients and the constant gravitational term, ρ 0 g, are zero), where f (x) is a function that expresses cos φ in terms of the coordinate x.Following Ref. [22], we assume that the projection of the wall law depends locally on the velocity, v, as , where G(v) has some generic properties defined below.
In order to close Eq.( 9) and Eq. ( 4) they need to be supplemented with proper constitutive equations for the temperature and solute concentrations.The main mechanisms that we consider for the evolution of the temperature are: • Convection (transport with the velocity), included in the second term of the left hand side of Eq. (4) • Heat exchange with an external source (Newton's cooling): Analogously, for the solute concentration, we consider the following physical mechanisms: • Convection (transport with the velocity), included in the second term of the left hand side of Eq. (5) • Solute diffusion: J S,1 = −c∇S • Onsager coupling (thermodiffusion or the Soret effect): J S,2 = b∇T .so from Eq. ( 5), Finally, we nondimensionalize the resulting equations: lengths are scaled by the loop length (so the whole loop is the interval [0, 1]), temperatures are scaled by (gα T ) −1 , solute concentration by (gα S ) −1 , . . .Finally, we arrive at our main system of equations The term l(v)(T a − T ) is the Newton's linear cooling law as in [16,15,14,28,17,29,5], which represents the heat transfer law across the loop, a positive quantity depending on the velocity, and T a is the (given) ambient temperature distribution, see [15,28,29,5].By physical consistency, in this paper we consider that the thermal diffusivity ν ≥ 0.
The system of equations ( 10) forms an ODE/PDE system for the evolution for the observables v, T and S. Finally, in the equation for the solute concentration S(t, x) the solute diffusivity, c > 0, and the Soret coefficient, b > 0, are both strictly positive, like in the previous models with this kind of binary fluids as [9,16,15,14,17].
Notice that, throughout this work, • dx = 1 0 • dx denotes integration along the closed path of the circuit.We can make this identification if we consider only periodic functions (with period 1).The function f describes the geometry of the loop and the distribution of gravitational forces [29,18], with f (x) dx = 0 (as f (x) = cos φ(x), so integrated along a closed loop it is zero).
The parameter ε in equation (10) is the nondimensional version of µ/E which has dimensions of time.Roughly speaking, ε gives the (nondimensional) time scale in which the transition from elastic to fluid-like occurs in the material.
We assume that G(v), which specifies the friction law at the inner wall of the loop, is positive and bounded away from zero.This function has been usually taken to be G(v) = G, a positive constant for the linear friction case [18] (Stokes flow), or G(v) = |v| for the quadratic law [11,20], or even a rather general function given by G(v) = g(Re)|v|, where Re is the Reynolds number, Re = ρvL/µ.Here we will consider a general function of the velocity assumed to be large [28,25].The functions G, f , and l incorporate relevant physical constants of the model, such as the cross sectional area, D, the length of the loop, L, the Prandtl, Rayleigh, or Reynolds numbers, etc., see [28].G and l as continuous functions, such that G(v) ≥ G 0 > 0, and l(v) ≥ l 0 > 0, for G 0 and l 0 are positive constants.Below, we will impose additional restrictions that will be introduced below in Eq. (28), this is we consider the nonlinear friction G such that G satisfies (28), i.e. there exists a constant h 0 ≥ 0 such that: 2.1.Well-posedness and boundedness: Global attractor.We will introduce some function spaces that will be used in the study of the existence of solutions of (10).Let Ω = (0, 1) and consider the spaces Ḣm per (0, 1) = H m loc (IR) ∩ L2 per (0, 1).( 13) 2.1.1.Existence and uniqueness of solutions.In this section, we prove the existence and uniqueness of solutions of the thermosyphon model (10), with f ∈ L2 per (0, 1), T a , T 0 ∈ Ḣ1 per (0, 1) and S 0 ∈ L2 per (0, 1), where L2 per (0, 1) and Ḣ1 per (0, 1) are given by ( 12) and we note that the dot stand for functions with zero average (and it is not related to time derivatives of the functions).For ν > 0, if we integrate the equation for the temperature along the loop, taking into account the periodicity of T , i.e., ∂T ∂x dx = ∂ 2 T ∂x 2 dx = 0, and Therefore, T (t, x)dx → T a (x)dx exponentially as the time goes to infinity for every T 0 (x)dx.
Moreover, if we consider τ (t, x) = T (t, x) − T (t, x)dx then from the second equation of the system (10), τ satisfies the equation: where τ a (x) = T a (x) − T a (x)dx.We integrate the equation for the solute concentration along the loop and taking into account the periodicity of S, ∂S ∂x dx = ∂ 2 S ∂x 2 dx = 0 and d dt ( S(t, x)dx) = 0.As S(t, x)dx is constant, it implies that the solute S(t, x) = S 0 (x) for all t.
Defining σ(t, x) = S(t, x) − S 0 (x)dx, (not be confounded with the stress component σ) then from the third equation of the system (10), σ satisfies the equation: Therefore, (v, τ, σ) satisfies the system (10) with τ a , τ 0 , σ 0 replacing T a , T 0 , S 0 respectively and f (x)dx = τ 0 (x)dx = τ a (x)dx = σ 0 (x)dx = 0 and T (t, x)dx = S(t, x)dx = 0 for all t ≥ 0. Hereafter we consider all the functions of the system (10) to have zero average.Also, as ν, c > 0 the operators νA = −ν ∂ 2 ∂x 2 and cA = −c ∂ 2 ∂x 2 , together with periodic boundary conditions, are unbounded, self-adjoint operators with compact resolvent in L 2 per (0, 1), that are positive when restricted to the space of zero average functions in L2 per (0, 1).Moreover, the equation for the temperature T and the equation for the solute concentration S in (10) are of parabolic type for ν, c > 0. We write the system (10) as the following evolution system for acceleration, velocity, temperature and solute concentration: That is: with and initial data ) and has compact resolvent, see Eq. (12).
Using the results and techniques about sectorial operator of [10] to prove the existence of solutions of the system, we have the Theorem 2.1.1.
Step (i).We prove the local existence and regularity.This follows easily from the variation of constants formula of [10].In order to prove this, we write the system as ( 15), and we have: where the operator per (0, 1) × Ḣ2 per (0, 1) and has compact resolvent.In this context, the operator A = − ∂ 2 ∂x 2 must be understood in the variational sense, i.e., for every T, ϕ ∈ Ḣ1 per (0, 1), and L2 per (0, 1) coincides with the fractional space of exponent 1 2 as in [10].We denote Ḣ−1 per (0, 1) as the dual space and .the norm on the space L2 per (0, 1).If we prove that the nonlinearity ) is well defined, Lipschitz and bounded on bounded sets, we obtain the local existence for the initial data in Y = IR 2 × Ḣ1 per (0, 1) × L2 per (0, 1).Using H(v) ≡ G(v)v and l(v) are locally Lipschitz together with f ∈ L2 per (0, 1) and T a ∈ Ḣ1 per (0, 1), we will prove the nonlinear terms, F 1...4 in Eqs. ( 16)-( 19) satisfy is well defined, Lipschitz and bounded on bounded sets.
Applying the techniques of variation of constants formula of [10], we obtain the unique local solution (w, v, T, S) ∈ C([0, t * ], Y) (with a suitable t * > 0) of ( 14), which are given by with where ) and using again the results of [10], (the smoothing effect of the equations together with the bootstrapping method), we get the regularity of solutions.
Step (ii).To prove the global existence, we must show that the solutions are bounded in Y = IR 2 × Ḣ1 per (0, 1) × L2 per (0, 1) for finite time intervals and using the nonlinearity of F, maps bounded on bounded sets, we conclude.
To prove that the norm of T is bounded in finite time, we multiply the equation for the temperature by T in L2 per (0, 1).Then integrating by parts, we have: 1 2 )dx = 0. Using Cauchy-Schwartz and Young inequalities and then the Poincaré inequality for functions of zero average, since T (t, x)dx = 0, and, as π 2 is the first nonzero eigenvalue of A = − ∂ 2 ∂x 2 in L2 per (0, 1), we obtain 1 2 and using the fact that l(v) ≥ l 0 > 0 we get and we conclude that the norm of T in L2 per (0, 1) remains bounded in finite time.By differentiating the third equation of ( 14) with respect to x, we obtain what proves that the norm of T in Ḣ1 per (0, 1) remains bounded in finite time.Then, we show that the norm of S in L2 per (0, 1) does not blow-up in finite time.Multiplying the fourth equation of ( 14) by S, integrating by parts, applying the Young inequality and again taking into account that S(t, x for every > 0 with C = 1 4 .Thus, taking = c 2 , and using (25) together with the Poincaré inequality for functions with zero average, we find with k 1 > 0. Therefore S(t) remains bounded in finite time.
Finally, since T and S are bounded in finite time, that implies that |w(t)|, |v(t)| remain also bounded in finite time.Hence we have a global solution in the nonlinear semigroup in Y = IR 2 × Ḣ1 per (0, 1) × L2 per (0, 1).2.1.2.Boundedness of the solutions: Global attractor.In this section, we adapt the results and techniques in Refs.[16,15,25] for a fluid with one component, to prove the existence of the global attractor for a binary fluid for the semigroup defined in the space Y = IR 2 × Ḣ1 per (0, 1) × L2 per (0, 1).To obtain the asymptotic bounds on the solutions as t → ∞, we consider the friction function G as in [16,15,25] satisfying the hypotheses of the previous section and there exits a constant h 0 ≥ 0 such that: Using the l'Hopital's lemma proved in [24] we have the following lemma proved in [32].
Lemma 2.1.2.If we assume that G(r) and H(r) ≡ rG(r) satisfy the hypotheses of Theorem 2.1.1,together with (28), then: with Remark 2.1.3.We note that the conditions (28) are satisfied for all the friction functions G considered in the previous works, i.e., the thermosyphon models where G is constant or linear or quadratic law.Moreover, the conditions (28) are true for G(s) ≈ A|s| n , as s → ∞.
Theorem 2.1.4.Under the above notation and hypotheses of Theorem 2.1.1,if we assume that G satisfies (29) for some constant H 0 ≥ 0 then Part (i) In particular: Part (ii) If ν = 0 and there exists L 0 a positive constant such that L 0 ≥ l(v) ≥ l 0 > 0, then for any solution of (10) in the space Y = IR 2 × Ḣ1 per (0, 1) × L2 per (0, 1) we have: ∂T a ∂x where ∂T a ∂x .
Finally, since the sectorial operator B, as defined in Theorem 2.1.1,has compact resolvent, the rest follows from [8] [Theorem 4.2.2 and 3.4.8].
Remark 2.1.5.First, we consider l(v) ≥ l 0 > 0 and we note that in Theorem 2.1.4,we assumed that also "l(v) ≤ L 0 ", this hypothesis is satisfied when we consider the Newton's linear cooling law l(v) = k(T a − T ), where k is a positive quantity i.e., l(v) = k = L 0 as [13].Moreover, this condition is also satisfied if we consider l = l(v)(T a − T ) where l(v) is a positive upper bounded function.
Second, it is important to note that we prove in the next section the existence of the global compact and connected attractor and the inertial manifold for the system (14), for the general Newton's linear cooling law, where l(v) is a locally Lipschitz function with l(v) ≥ l 0 > 0 without the additional above hypothesis for l(v), namely, "without the upper bounded L 0 ", and for every ν ≥ 0.
In order to get this, we introduce the Fourier expansions of the fields T and S to prove lim sup t→∞ T (t) ≤ T a and lim sup t→∞ S(t) ≤ b c T a for every locally Lipschitz and positive function l(v) and for every ν ≥ 0.
2.2.Asymptotic behavior: Reduction to finite-dimensional systems.In this section we reduce the asymptotic behavior of the system (14) to study the other finite dimensional systems in some cases.We take a close look at the dynamics of ( 14) by considering the Fourier expansions of all the functions of the system (temperature and solute concentration).First, we consider the Fourier expansion for the function associated to the geometry of the loop f, and the ambient temperature T a , whose coefficients are very important to study the asymptotic behavior of the system.Then, we prove the asymptotic behavior of the system (14), described by suitable Fourier coefficients associated to f and T a .
We note the Fourier expansion for all g ∈ Ḣm per (0, 1), m ≥ 0 is given by the expression g(x) = k∈Z * a k e 2πkix with Z * = Z \ {0} and we have Assume that T a , T ∈ Ḣ1 per (0, 1) and f, S ∈ L2 per (0, 1) are given by the following Fourier series expansions: with the initial data T 0 ∈ Ḣ1 per (0, 1) is given by T 0 (x) = k∈Z * a k0 e 2πkix and S 0 ∈ L2 per (0, 1) is given by S 0 (x) = k∈Z * d k0 e 2πkix .Since all functions involved are real and periodic, we have āk Proposition 2.2.1.Under the above notation and hypotheses of Theorem 2.1.1,we consider T a ∈ Ḣ1 per (0, 1) and f ∈ L2 per (0, 1) given by (53) and the initial data T 0 ∈ Ḣ1 per (0, 1) given by T 0 (x) = k∈Z * a k0 e 2πkix and S 0 ∈ L2 per (0, 1) given by S 0 (x) = k∈Z * d k0 e 2πkix .Let (w, v, T, S) be the solution of the system (10) given by Theorem 2.1.1;then we have: (i) The coefficients a k (t) and d k (t) in (54), satisfy the equations The equation for the velocity is Proof.By using the following Fourier expansion for a generic function u(t, x) of the form, u(t, x) = k∈Z * u k (t)e 2πkix with Z * = Z \ {0} then, the Fourier coefficients u k (t) satisfy the following relations: Consider the model (10) and the Fourier series expansions of all functions depending on spatial variable x, i.e.
given by (54), T a (x) = k∈Z * b k e 2πkix and f (x) = k∈Z * c k e 2πkix given by (53), together with the expansions of initial data for temperature T 0 (x) = k∈Z * a k0 e 2πkix and for solute S 0 (x) = k∈Z * d k0 e 2πkix ; then we easily find that the coefficients for temperature a k (t) and solute concentration d k (t) are the solution of (55).Besides, it is sufficient to note that Remark 2.2.2.We note that the system (10) is equivalent to the system (14) for acceleration, velocity, temperature and solute concentration and from the above proposition, it is equivalent to the following infinite system of ODEs (60) The system of equations (60) reflects two of the main features: (i) the coupling between the modes enter through the velocity, while the diffusion acts as a linear damping term, (ii) the non linear term, given by Newton's cooling law in the temperature also affects the evolution of the solute concentration.
In what follows, we will exploit this explicit equation for the temperature and solute concentration Fourier modes to analyze the asymptotic behavior of the system and to obtain the explicit low-dimensional models.A similar explicit construction was given by Bloch and Titi in [2] for a nonlinear beam equation where the nonlinearity occurs only through the appearance of the L 2 norm of the unknown.A related construction was given by Stuart in [26] for a nonlocal reaction-diffusion equation.
In the following section, we prove the boundedness of these coefficients which improve, in some sense, the boundedness of temperature and solute concentration, in order to prove the existence of the inertial manifold for the system (14), using similar techniques as in Refs.[22,25].
2.3.Inertial manifold.We consider the general case ν > 0 together with the nonlinear Newton's cooling law l(v)(T a − T ) with l(v) ≥ l 0 > 0 being locally Lipschitz function and use inertial manifold techniques, in the spirit of the non-diffusion case of [22], to give an explicit low-dimensional system of ODEs that describes the asymptotic dynamics of (14).The existence of an inertial manifold does not rely, in this case, on the existence of large gaps in the spectrum of the elliptic operator but on the invariance of certain sets of Fourier modes.Proposition 2.3.1.Under the above notation and hypotheses of Theorem 2.1.4,with initial conditions (w 0 , v 0 , T 0 , S 0 ) ∈ Y = IR 2 × Ḣ1 per (0, 1) × L2 per (0, 1), for every solution of the system ( 14), (w, v, T, S), and for every k ∈ Z * , and recalling the expansions with the initial data T 0 ∈ Ḣ1 per (0, 1) is given by T 0 (x) = k∈Z * a k0 e 2πkix and S 0 ∈ L2 per (0, 1) is given by S 0 (x) = k∈Z * d k0 e 2πkix , we have In particular, if T a ∈ Ḣm per (0, 1), we have the global compact and connected at- Proof.(i) From (55), we have (69) and taking into account that we obtain: From (55), we have Substituting (69) in (72), we have where ]dr e −4νπ 2 k 2 (s−r) dr ds.
Using (70) and l(v) ≥ l 0 in I 1 (t), we get, Again, as |e ))e − s r l(v) dr ds so using Eq. ( 70), we find Finally, from (73) together with (74) and (75), we have (ii) Using again Theorem 2.1.4and taking into account Eq. ( 52) together with we obtain for any solution of ( 14), we have (65), (66), ( 67) and (68), since the sectorial operator B define above (in the section 2.1.1)has compact resolvent from [Theorem 4.2.2 and 3.4.8] in Ref. [8], the system has a global compact and connected attractor,A, in 61) and (62), for any (w(t), v(t), T (t, x), S(t, x)) ∈ A we have per (0, 1) then T, S ∈ Ḣm per (0, 1) and we have (w(t), v(t), T (t, x), S(t, x)) Finally we show that C is compact in Ḣm per (0, 1).Indeed, for any sequence {R n } in C we can extract a subsequence that we suntil denote {R n } such that it converges weakly to a function R and such that, for any k ∈ Z * , the Fourier coefficients verify r n k → r k as n → ∞, where r k is the kth Fourier coefficients of R. Therefore, |r k | ≤ d|b k | and for every integer N 0 .
where .m denotes the norm in Ḣm per (0, 1).Hence, the first term goes to zero as n → ∞ and the second can be made arbitrarily small as ) and the result is proved.Now, we will prove that there exists an inertial manifold M (see a definition in Ref. [7]) for the semigroup S * (t) in the phase space Y = R 2 × Ḣ1 per (0, 1)× L2 per (0, 1), i.e., a submanifold of Y such that (i)S * (t)M ⊂ M for every t ≥ 0, (ii) there exists δ > 0 satisfying that for every bounded set B ⊂ Y, there exists C(B) ≥ 0 such that dist(S(t), M) ≤ C(B)e −δt , t ≥ 0 see, for example, [7] and [23].
Assume that T a ∈ Ḣ1 per (0, 1) with with b k = 0 for every k ∈ K ⊂ Z * with 0 / ∈ K, since T a (x)dx = 0. We denote by V 1 and V 0 the closure of the subspaces of Ḣ1 per (0, 1) and L2 per (0, 1) respectively generated by {e 2πkix , k ∈ K}.
From (71), and taking into account that b k = 0 for k ∈ Z * , we have |a k (t)| ≤ |a k0 |e −4νπ 2 k 2 t and together with 4νπ 2 k 2 t ≥ 4νπ 2 t for every k ∈ Z * with (52), implies that T 2 (t) . Then, using (α + β + γ) 2 ≤ 4(α 2 + β 2 + γ 2 ), together with π 2 k 2 δt ≥ π 2 δt for every k ∈ Z * and δ = min{ν, c}, we get From this, we obtain Now, we note for any > 0 there exists N 0 ( ) such that |I k − 1| ≤ for every |k| ≥ N 0 , and then |I k | 2 ≤ ( + 1) 2 .Thus taking = 1 (for example), we have with c k = 0 for every k ∈ J ⊂ Z and c k = 0 if k / ∈ J.. On the inertial manifold So, the evolution of velocity v, and acceleration w depends on the Fourier coefficients of T and S in the set K ∩ J.In order to solve full system, this, first we get the coefficients of T and S which belong to the set K ∩ J and after we must solve the equations for the coefficients of T and S k / ∈ K ∩ J (i.e.|K \ (K ∩ J)|).We note that 0 / ∈ K ∩ J and since K = −K and J = −J then the set K ∩ J has an even number of elements, that we denote by 2n 0 .Therefore, the number of positive elements of K ∩ J, (K ∩ J) + , is n 0 .
Proof.On the inertial manifold Therefore, the dynamics of the system depends on the coefficients in K ∩J.Moreover the equations for a −k and d −k are conjugates of the equations for a k and d k , and therefore we have From this, taking real and imaginary parts of a k , (a k 1 , then on the inertial manifold we get a homogeneous linear equation for the velocity with positive coefficients, that is: and therefore we have v(t) → 0 when t → ∞.Moreover from the equation for the temperature in (14) we have that the function θ = T − θ ∞ satisfies the equation: We can multiply by θ in L2 per and taking into account that ∂θ ∂x θdx = 1 2 and using Cauchy-Schwartz and Young inequality with δ, C δ = 1 4δ and then the Poincaré inequality, since θdx = 0, we have 1 2 Using v(t) → 0, together with Gronwall lemma, we prove that θ(t) → 0 in L2 per (0, 1).We note that S − b c θ ∞ satisfies the equation multiplying this equation by S − b c θ integrating by parts and taking into account, the periodicity of S, we obtain ∂S ∂x Sdx = 0, and applying the Young and Poincaré inequality, we get the inequality du/dt + c * u ≤ 2 for every t ≥ t 0 , large enough, with u = S − b c θ ∞ 2 , and c * > 0. Finally, from the Gronwall lemma we have and thus lim sup t →∞ u(t) ≤ 2 for every > 0.
Remark 2.4.2.Taking the real and imaginary parts of the coefficients of the temperature, a k (t), heat flux at the wall of the loop, b k , the geometry of the circuit, c k , and the solute concentration, the asymptotic behavior of the system ( 14) is given by a reduced explicit system of ODEs in IR N , with N = 4n 0 + 2, given by + Thus, we have reduced the asymptotic behavior of the initial system (14) to the dynamics of the reduced explicit system (87).It is worth noting that, from the above analysis, it is possible to design the geometry of the circuit and/or the external heating source, by properly choosing the functions f and/or the ambient temperature, T a , so that the resulting system has an arbitrary number of equations of the form N = 4n + 2.
Note that K and J may be infinite sets, but their intersection is finite.For instance, for a circular circuit we have f (x) ∼ a sin(x) + b cos(x), i.e., J = {±1} and then K ∩ J is either {±1} or the empty set.
3. Numerical integration.In this section we present some numerical integrations using the MATHEMATICA package [31] for the resolution of the differential equations (87), using a fourth-order explicit Runge-Kutta method suitable for stiff equations [5,10].To study the asymptotic behavior of the system of equations, we consider the coefficients of temperature a k (t) and the coefficients of solute concentration d k (t) with k ∈ K ∩ J, as the relevant modes.Then, we have the finite system of differential equations: (87) If we take, for sake of simplicity, a circular geometry, then J = {±1} and K ∩J = {±1}.Also, if we take k = 1 and omit the equation for −k, the conjugate of k.Therefore, we have the following transformed set of equations: where the unknowns are w(t), the acceleration of the fluid, v(t), the velocity of the fluid, a 1 (t), the Fourier mode of the temperature and d 1 (t), the Fourier mode of the solute concentration.
In order to reduce the number of free parameters, we make the following change of variables a 1 c −1 → a 1 and d 1 c −1 → d 1 and we introduce the real and imaginary parts of the equations in the following way: This forms a system of six differential equations with six unknowns where we need to make explicit choices for the constitutive laws of both the fluid-mechanical and thermal properties.The function G(v) has a clear physical meaning; it interpolates between a low Reynolds number friction law (in which the overall friction, G(v)v is nonlinear (Stokes friction law), and high Reynolds number (in which the friction is a quadratic law).For the friction law G(v) and heat flux l(v) = l(v)(T a − T ), we will take the ones used in the references [5,10].For the numerical experiments, which are of a similar model of thermosyphon for a fluid with one component, they use the functions G(v) = (|v| + 10 −4 ) and l(v) = (10 −2 |v| + 1).
The parameter ε in the system of equations ( 14) is the nondimensional version of µ/E which has dimensions of time.Roughly speaking, it gives the (nondimensional) time scale in which the transition from elastic to fluid-like occurs in the fluid.Numerical analysis has been carried out keeping ε the viscoelastic coefficient as the tuning parameter ranging from 10 −5 (almost Newtonian) to 10 1 a highly viscoelastic value.The impact of ε on the system has been keenly observed for the parameters of time t, as short as 100 time units and as long as 1000 time units.All the variables and equations that we deal with are non-dimensional.
For the Soret effect diffusion coefficients (b and c), we will assume the values calculated by Hart in [9], that consider a thermosyphon of circular geometry of radius R 0 (for the loop) and R p (for the pipe).Hart takes the values for a mixture of alcohol and water, borrowed from Hurle and Jakeman [9].This reference settles down that c = Ds V R0 is the number of Lewis, where D s is the diffusivity of the solute that has a value for such a mixture of 10 −5 cm 2 s −1 and V is the scale of velocity, with a value of 10 −2 cms −1 for a circular thermosyphon whose loop to pipe radius ratio is 10.Therefore we will take c = 10 −3 .Also, as Hart discussed in [9], b (the Soret diffusion coefficient) is a parameter that determines the qualitative behaviour of the variables.Therefore, in the numerical experiments we treat the value of b the Soret coefficient as another tuning parameter ranging from 10 −5 to 10 1 .
We plot the relevant modes of temperature and solute concentration as they are the ones to have influence on the acceleration w(t) and velocity v(t).The real and imaginary parts of the Fourier transform of the temperature (a 1 and a 2 ) and the analogous for the solute concentration (respectively, d 1 and d 2 ).As the system is multidimensional, we present the results in temporal graphs (a given variable versus time) and phase-space graphs (two physical variables plot against each other).We will show that, in analogy with the classical Lorenz system, as ε varies, the dynamics of the model undergoes various transformations including steady asymptotic behavior, meta-stable chaos, i.e., transient irregular behavior followed by convergence to equilibria, periodic behaviors and chaotic progressions.4. Results and discussion.In this section we analyze and discuss the behavior of the system for various parameters of the viscoelastic coefficients with various Soret Effect gradients.In each of the subsequent subsections, we define and fix the various parameters that are employed in this model.A and B refer to the positiondependant (x) heat flux inside the loop.Without loss of generality, we will assume A = 0, while B is fixed to be 30, in order to simplify and, in analogy with Lorenz's model, as it is shown in references [5,10] (changing A and B simultaneously only results in a change in the phase of initial temperature profile).The heat diffusivity is fixed to ν = 0.002, namely, temperature diffusion is present but not dominant.The initial conditions are fixed as w(0) = 0, v(0) = 0, a 1 (0) = 1, a 2 (0) = 1, d 1 (0) = 0.01, and d 2 (0) = 1.As, in all cases, we are concerned with the long time behavior of the system, the specific choice of the initial condition does not change significantly the main qualitative conclusions following our results.Clearly, the chaotic behavior of the system is observed for these ranges of parameters.The plots in Fig. 2 of acceleration and velocity portray the chaotic nature of the behavior of the system.The acceleration ranges from -20 to 20 while the velocity ranges from -3 to 3. As the velocity is the time derivative of acceleration, the deviation in the time series plot of velocity is reduced to -3 to 3.But both the plots exhibit chaotic progresses for the entire range of 1000 time units.So, when the viscoelastic coefficient as well as the Soret coefficient are small enough (∼ 0.00001) the system has, apparently, chaotic behaviors.
The real and imaginary parts of the temperature are plotted in Fig. 3 for the last time units used in the numerical integration.Both plots exhibit complex and chaotic  features until the end of the entire 1000 time units.As in the case of acceleration and velocity, the temperature too is chaotic, when both the viscoelastic and Soret effect coefficients are as small as 10 −5 .This is a consequence of the coupling between the variables.The right panel in Fig. 3 shows the phase-space graph, resembling the celebrated Lorenz's phase-diagram.
The real and imaginary parts of the solute concentration Fourier mode are given in Fig. 4 for the last time units, the first one corresponding to the real part of the solute concentration and, the second, to the imaginary part of the solute concentration.Both plots have a similar chaotic behavior.Comparing the plots of temperature with solute concentration, we find that the deviation in temperature is very well distributed throughout the progression, whereas the deviation in the solute concentration is found to be more concentrated at the centre.
To sum up this first numerical experiment, we find that when the value of viscoelastic and Soret coefficients are both 10 −5 , the system is chaotic.Similar type of chaotic behaviors are found for all the values of the viscoelastic components ε ranging from 10 −5 to 10, when the value of Soret coefficient is 10 −5 .From the above observations, we can qualitatively state that when the Soret coefficient b is smaller than 10 −5 , the system has chaotic behaviors.In Sec.4.4 we will extend these results quantitatively by means of standard Lyapunov exponent analysis.The phase-space plots of system variables is shown in Fig. 5, portray the transition that takes place in the nature of the behavior of the system, though it suntil is chaotic apparently.The plot exhibits a chaotic behavior at the early stages of the dynamics together with a decay towards the long-time behavior.Thus, after around 400 time units, the system begins to change its behavior to converge to a periodic solution (not shown).This change of behavior is present also in the time series plots of temperature and solute concentration, which help us to make note of the changing phenomena for these particular parameters.4.2.2.Case II.The quasi-periodic behavior for ε=1.Continuing from the previous numerical experiment, when ε is increased to 1, a relatively greater value from the previous experiment, we find a significant change to quasi-periodic behavior taking place in the system.In order to illustrate the quasi-periodic behavior of the system, we plot different phase-diagrams in Fig. 6.In contrast with the phase-space plots in Fig. (4) in this case the solution covers the ring in that figure, always in the same direction (either clockwise or anti-clockwise) what shows that the behavior is (quasi-)periodic but it is no longer chaotic.
To sum up, we have found that when the viscoelastic coefficient is 1 and the Soret coefficient b = 10 −3 , the system displays a transition from a chaotic behavior to a quasi-periodic behavior.It also shows that the system is slowly transforming towards a state of stable behavior.4.2.3.Case III.The stable behavior for ε=10.Finally, numerical experiments were carried out for the viscoelastic coefficient ε=10, keeping the Soret coefficient fixed to be b = 10 −3 .For ε=10, a relatively large value, the velocity displays a chaotic beginning but as time progressed it begins to stabilize around the time unit 400 (not shown).Once crossing the 400 time mark, the plot shows a very stable progress until the end of 1000 time units.A transition from chaotic to stable behaviors of the system is displayed by the variables a 1 , a 2 , d 1 and d 2 (not shown).This gives us a broad understanding that when the viscoelastic coefficient is large i.e., ε=10, the system tends to stabilize as the time progresses.This result can be confirmed by the value of the maximum Lyapunov exponent for this set of parameters and summarized in Table 2.In particular, for ε = 10 and b = 10 −3 , this exponent is −0.03 so, as it is the maximum, all the Lyapunov exponents are negative and the system converges to a stable solution.
To sum up, we find that when the value of viscoelastic component is 10 and the Soret coefficient is 0.001, the system has a transformation from chaos to stability.From the above observation, we can state that when the value of viscoelastic component is 10, the system has stabilizing effects.
To sum up the numerical experiment II, we find that when the Soret coefficient b = 10 −3 , the system has different kind of behaviors depending on the range of ε the viscoelastic coefficient.When ε is relatively small, i.e., ε < 10 −1 the system has chaotic effect.But when ε=1, the system has quasi-periodic effect and at ε=10, it begins to stabilize.From these results, we can state that when the Soret coefficient b = 10 −3 , the system has different kind of behaviors depending on the values of ε the viscoelastic coefficient.4.3.Experiment III: Soret coefficient b = 1.In the third set of numerical experiments, numerical integrations were carried out, keeping the Soret coefficient b fixed to be 1, while the viscoelastic coefficient ε ranging from 10 −5 to 10 as tuning parameter.A significant change in the behavior of the system is observed for the viscoelastic coefficient ε=1.That is, when the viscoelastic coefficient and the Soret coefficient are equal to 1, the system has a stable behavior.
In all cases integrated, we can claim that, in general, the most important parameter is ε.
Concluding all the three sets of numerical experiments, we want to emphasize that the overall impact of the viscoelastic coefficient on the system for the entire range of parameter ranging from 10 −5 to 10  1. Qualitative summary of the overall behavior of the system for different values of the viscoelastic characteristic time, ε (columns) and the Soret gradient(rows).We introduce the following notation to account for the obtained numerical results: 'S' a stable behavior,'C' denotes a fully chaotic behavior, and 'QP' a transitional behavior from chaotic to either periodic or quasiperiodic.
Soret gradients.We observe that when the viscoelastic coefficient is small, i.e., ε = 10 −5 , 10 −4 , 10 −3 and 10 −2 the system has chaotic behaviors irrespective of the Soret gradient.But when ε=0.1, 1 and 10 the system has three different kind of behaviors depending on the Soret gradients.At first, when the Soret coefficient is 10 −5 and 10 −4 , the system has chaotic behavior for all the values of viscoelastic coefficients.In the second place, when the Soret gradient is 10 −3 , the system has chaotic, quasi-periodic and stable behavior when ε=0.1, 1 and 10 respectively.And in the third place, when the values of viscoelastic coefficient and the Soret gradient are 0.1, 1 and 10, the system has stable behavior.Table 1 provides the qualitative details of all of these behavior of the system for various parameters, where we have introduced the following notation to account for the obtained numerical results: 'S' a stable behavior,'C' denoting a fully chaotic behavior, and 'QP' a transitional behavior from chaotic to quasi-periodic behavior.4.4.Analysis of the behavior of the system using Lyapunov exponents.The behavior of physical systems has been a matter of primary importance to scientists, engineers and mathematicians in order to determine and characterize the dynamical behavior of the system.In this regard, the behavior of the system that we are working with has to be ascertained as it has nonlinear dynamics with many variables.To that aim, the Russian mathematician Alexandr Lyapunov introduced definitions and criteria to establishes unambiguously chaotic, periodic, quasi-periodic or stable behavior by studying the linearization of the equations of motion to determine the behavior of any system, through the now so-called Lyapunov exponents.The signs and the values of the Lyapunov exponents allow us to determine the qualitative and quantitative patterns of behavior of any system [30].
The Lyapunov exponents have been proved useful for determining and distinguishing the various types of orbits and behavior of our system.In the first case, when the viscoelastic coefficients ε = 10 −5 , 10 −4 , 10 −3 and 10 −2 relatively lesser gradients, the system has positive Lyapunov exponents, which is a clear indicator for chaotic or strange behavior of the system.In the second case, when the viscoelastic coefficients ε are 0.1, 1 and 10, relatively higher gradients, the system has positive, negative and around the zero point Lyapunov exponents which imply different kind of behaviors.In this case, depending on the value of Soret gradients the system exhibits chaotic, quasi-periodic, periodic and stable behaviors.As the viscoelastic and Soret gradients increased, the system has more negative Lyapunov exponents which confirm to greater stable behaviors.In a particular case, when the viscoelastic gradient is 1 and the Soret gradient is 10 −3 , the system has the Lyapunov exponent in and around zero, which is an indication for a quasi-periodic behavior.2. The maximum Lyapunov exponent of the system for different values of the viscoelastic characteristic time, ε (columns) and the Soret gradient(rows).We can assume that maximum Lyapunov exponents close to 0±0.1 correspond to quasi-periodic behavior (as simple inspection of the time series plots confirm).As Table 2 shows, we can observe that for a major part of the experiment, chaos is the common feature of the behavior of the system, apart from a few stable and quasi-periodic behavior of the system where the viscoelastic as well as the Soret coefficient are high.Figure 7 gives the overall behavior of the system, where the dark area indicates the chaos, the shaded area indicates the stable behavior and the white area indicates the quasi-periodic behavior of the system.

Conclusion.
In this work we have derived a system of equations that governs the motion of a viscoelastic material with the Soret effect inside a closed-loop thermosyphon.This model, which is a generalization of a previous work [32], has been used to study the presence of complex chaotic behaviors and also to relate them with the underlying viscoelastic (memory effects).This generalization comes from the physical (as it allows to study the coupling of new physical mechanisms, as the Soret effect) and mathematical sides (as the new field describing the solute concentration is coupled in a non-trivial way with the velocity of the fluid and, also, with the temperature).Thus, new equations for the solute concentration, S, introduces new complexity into the equations and, consequently, the results for the inertial manifold cannot be straightforwardly derived from Ref. [32] but, rather, involve new assumptions and calculations with the Fourier coefficients as shown in Sec. 2.
The results suggest that, when the value of the viscoelastic coefficient ε is small, it drives the dynamics to a chaotic behavior for all the physical observable (acceleration, velocity, temperature and solute concentration).In contrast, and contrary to what happens in the absence of solute, as the value of ε gradually increases, the system has a transition from chaotic to quasi-periodic or stable behavior depending on the Soret gradient.From these numerical experiments, we observe that the viscoelastic material with Soret gradients sets the system in chaos when the viscoelastic coefficient ε is small and stabilizes the system when the viscoelastic coefficient ε is higher.
Physically, this means that, when the viscoelastic effects are large (namely, the time scale ε −1 is comparable with the characteristic time scale of variation due to thermodiffusion and Soret effect, the memory smoothening arising from equation (3) drives the system towards an stable fixed point.This might explain why chaotic behaviors are more commonly observed in fluids than in solids in which (sustained) elastic oscillations are periodic or damped out with time due to dissipation.
This result prevents the system to develop chaotic oscillations, what provides a good candidate for engineering tuning of the coolant properties of the fluid inside the thermosyphon.Thus, the inclusion of a highly thermophoretic solute allows the oscillations to damp-out and the system can sustain a constant (stable) velocity, thus increasing the performance of the device.
since all the functions involved are real and periodic, we have for allk ∈ Z * = Z\{0}, āk = a −k , bk = b −k , ck = c −k and dk = d −k ; to conclude that (10) is equivalent to infinite systems of ODEs consisting of (55) coupled with where M, N are the upper bounds for acceleration and velocity as given in (68) and (67) respectively and T, S ∈ C = {R(x) ∈ Ḣm per (0, 1), R(x) = k∈Z * r k e 2πkix , |r k | ≤ d|b k |}, where d = max{1, b c } and A is a compact set in IR 2 × Ḣm per (0, 1) × Ḣm per (0, 1).

4. 1 .
Experiment I: Soret coefficient b = 10 −5 .In this very first numerical experiment, keeping the Soret coefficient b as small as 10 −5 , we observe the impact of ε the viscoelastic coefficient on the system.Numerical experiments were carried out for different viscoelastic coefficients ε ranging from 10 −5 to 10 as tuning parameters.In Figs.1-3we show the behavior of the system for a particular case of numerical experiments for the values of ε = 10 −5 and Soret coefficient 10 −5 , plotting acceleration, velocity, temperature and the solute concentration versus time.

4. 2 .
Experiment II: Soret coefficient b = 10 −3 .In this second set of numerical experiments, by keeping the Soret coefficient b = 10 −3 , we determine the impact of ε the viscoelastic coefficient on the system, by ranging from ε = 10 −5 to 10.

4. 2 . 1 .
Case I. Transition from chaotic behavior for ε = 10 −1 .Numerical experiments were carried out for the viscoelastic coefficients 10 −5 , 10 −4 , 10 −3 , 10 −2 , and 10 −1 , keeping the Soret coefficient fixed to be b = 10 −3 .Until the value of ε=0.1, the system has chaotic behavior.But at ε=0.1, the system begins to change from fully chaotic behavior to periodic or stable behavior.To avoid redundancy, we will show, hereafter, the plots for the velocity.

Figure 7 .
Figure 7.The overall behavior of the system for different values of viscoelastic and Soret coefficients (the red area at the bottom indicates chaos, the blue area in the middle indicates the quasiperiodic behavior and the yellow area at the top indicates the stable behavior).