Abstract

Fluids subject to thermal gradients produce complex behaviors that arise from the competition with gravitational effects. Although such sort of systems have been widely studied in the literature for simple (Newtonian) fluids, the behavior of viscoelastic fluids has not been explored thus far. We present a theoretical study of the dynamics of a Maxwell viscoelastic fluid in a closed-loop thermosyphon. This sort of fluid presents elastic-like behavior and memory effects. We study the asymptotic properties of the fluid inside the thermosyphon and the exact equations of motion in the inertial manifold that characterizes the asymptotic behavior. We derive, for the first time, the mathematical derivations of the motion of a viscoelastic fluid in the interior of a closed-loop thermosyphon under the effects of natural convection and a given external temperature gradient.

1. Introduction

Instabilities and chaos in fluids subject to temperature gradients have been the subject of intense work for its applications in engineering and in atmospheric sciences. In such sort of systems, the fluid displays nontrivial behaviors (as turbulence or the formation of convective rolls) when subject to a heating that competes with buoyancy effects. A traditional approach that goes back to the pioneering work by Lorenz consists of the study of the system under some simplifications. Another approach is to study the controlled setups that capture the underlying complexity of the full system, being a thermosyphon one of those simpler cases [1].

In the engineering literature, a thermosyphon is a device composed of a closed-loop pipe containing a fluid where some soluble substance has been dissolved [2, 3]. The motion of the fluid is driven by the action of several forces such as gravity and natural convection. The flow inside the loop is driven by an energetic balance between thermal energy and mechanical energy. The interest on this system comes both from engineering and as a toy model of natural convection (for instance, to understand the origin of chaos in atmospheric systems).

Thus far, all the works on thermosyphons analyze the behavior of a Newtonian fluid inside the loop, consequently neglecting elastic effects in the system coming from either the fluid itself or the elastic walls of the loop. However, many interesting fluids are known to behave slightly different from the common (Newtonian) fluids in terms of their response to an applied stress and are commonly referred to as viscoelastic. Among them, it is worth emphasizing volcanic lavas, snow avalanches, flowing paint, or biological mucosas membranes.

Here, we consider a thermosyphon model in which the confined fluid is viscoelastic. This has some a priori interesting peculiarities that could affect the dynamics with respect to the case of a Newtonian fluid. On the one hand, the dynamics has memory and so its behavior depends on the whole past history, and, on the other hand, at small perturbations the fluid behaves like an elastic solid and a characteristic resonance frequency could, eventually, be relevant (e.g., consider the behavior of jelly or toothpaste).

The simplest approach to viscoelasticity comes from the so-called Maxwell constitutive equation [4]. 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:

For a Maxwellian fluid, the stress tensor takes the form where is the fluid viscosity, the Young’s modulus, and the shear strain rate (or rate at which the fluid deforms). Under stationary flow, (2) reduces to Newton’s law, and consequently, (1) reduces to the celebrated Navier-Stokes equation. On the contrary, for short times (where impulsive behavior from rest can be expected) equation (2) reduces to Hooke’s law of elasticity.

Memory effects can be well understood from (2) after performing a separation of variables and integrating. Thus, we can rewrite where it is clear that the local state of stress is calculated from the present and values of with a memory time scale of order .

In a thermosyphon, the equations of motion can be greatly simplified because of the quasi-one-dimensional geometry of the loop. Thus, we assume that the section of the loop is constant and small compared with the dimensions of the physical device, so that the arc length coordinate along the loop () gives the position in the circuit. The velocity of the fluid is assumed to be independent of the position in the circuit; that is, it is assumed to be a scalar quantity depending only on time. This approximation comes from the fact that the fluid is assumed to be incompressible, and so besides the quasi-one-dimensional assumption. On the contrary, temperature is assumed to depend both on time and position along the loop.

The derivation of the thermosyphon equations of motion is similar to that in [2, 3]. The simplest way to incorporate (2) is by differentiating (1) with respect to time and replacing the resulting time derivative of with (2). This way to incorporate the constitutive equation allows to reduce the number of unknowns (we remove from the system of equations) at the cost of increasing the order of the time derivatives to second order.

The resulting second-order equation is then averaged along the loop section (as in [2]). Hence, after nondimensionalizing the variables (to reduce the number of free parameters) we arrive at our main system of equations where ; that is, we consider Newton’s linear cooling law as in [3, 510] also in this paper we consider the diffusion of temperature given by the term . The parameter in (5) is the nondimensional version of 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.

The model we will take into consideration forms an ODE/PDE system for the velocity , the distribution of the temperature of the fluid into the loop, (5) with , where denotes integration along the closed path of the circuit. The function describes the geometry of the loop and the distribution of gravitational forces [2, 3]. Note that .

The system of (5) is not a trivial extension of the Newtonian model in [11] due to the first term in the differential equation for the velocity. Specifically, the addition of a term proportional to the second derivative of is singular, in the sense that it changes qualitatively the character of the equations. The implications of this singular perturbation cannot be ascertained using a standard boundary layer analysis (see the Appendix for details). So a complementary theoretical and numerical approach is mandatory.

We assume that 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 , a positive constant for the linear friction case [2] (Stokes flow), or for the quadratic law [12, 13], or even a rather general function given by , where is the Reynolds number, . Here we will consider a general function of the velocity assumed to be large [9, 14]. The functions , , and incorporate relevant physical constants of the model, such as the cross-sectional area, , the length of the loop, , and the Prandtl, Rayleigh or Reynolds number; see [9]. Here, we consider Newton’s linear cooling law where represents the heat transfer law across the loop wall, and is a positive quantity depending on the velocity, and is the (given) ambient temperature distribution, see [3, 5, 9, 10].

Hereafter, we consider and as continuous functions, such that , and , for and positive constants.

Our contributions in this paper are the following.(i)To obtain the system of (5) governing a closed loop thermosyphon model with a viscoelastic fluid and to study the asymptotic behavior of this system which is a generalization of the previous model [11]. Thus, although the dynamics of viscoelastic fluids has been extensively studied in the engineering literature, it has scarcely been considered in the applied mathematics community.(ii)To present an analysis beginning with the well-posedness and boundedness of solution. The existence of an attractor and an inertial manifold is shown and an explicit reduction to low-dimensional systems is obtained. It is noteworthy that we are able to obtain an exact finite-dimensional reduction (87) that may have a much lower number of degrees of freedom.(iii)To provide a detailed numerical analysis of the behavior of acceleration, velocity, and temperature which includes a thorough study of the various behaviors of the system for different values of viscoelastic fluid and ambient temperature distribution.(iv)The numerical analysis will show that viscoelasticity induces a chaotic behavior that is not captured by a boundary layer analysis (that would predict the same qualitative behaviors as in the original model in [11] see the Appendix for details) being the new (nontrivial) emergent behaviors induced by the viscoelasticity worth characterizing. The structure of the paper is as follows. The first section provides a general introduction to the system, explaining briefly the dynamics of a thermosyphon, viscoelastic fluids, and the objectives of this work. In Section 2, we prove the existence, uniqueness, and boundedness of the solutions. In Section 3, we provide a detailed derivation of the dynamics of the system in the inertial manifold as a reduced dimensionality version of the full system of (5). In Section 4, we present the numerical integration of the reduced system of equations valid in the manifold to understand the role of the main parameters of the physical system. Finally, in Section 5, we conclude with a proposal for future works.

2. Well-Posedness and Boundedness: Global Attractor

2.1. Existence and Uniqueness of Solutions

In this section we prove the existence and uniqueness of solutions of the thermosyphon model (5).

First, we observe that for , if we integrate the equation for the temperature along the loop taking into account the periodicity of , that is, , we have . Therefore, exponentially as time goes to infinity for every .

Moreover, if we consider , then from the second equation of system (5), we obtain that verifies the equation where .

Finally, since , we have and the equation for reads

Thus, with , we get verifying the system (5) with replacing , respectively, and now . Therefore, hereafter we consider the system (5) where all functions have zero average.

Also, if , the operator , together with periodic boundary conditions, is an unbounded, self-adjoint operator with compact resolvent in that is positive when restricted to the space of zero average functions . Hence, the equation for the temperature in (5) is of parabolic type for .

2.1.1. The Case with Diffusion:

We consider the acceleration and write the system (5) as the following evolution system for the acceleration, velocity, and temperature: that is, with and and the initial data .

The operator is a sectorial operator in with domain and has compact resolvent, where

Using the results and techniques of sectorial operator of [15], we obtain Theorem 1.

Theorem 1. We assume that and are locally Lipschitz, , , and . Then, given , there exists a unique solution of (5) satisfying where and for every . In particular, (5) defines a nonlinear semigroup, in , with .

Proof. We cover several steps.
Step 1. We prove the local existence and regularity. This follows easily from the variation of constants formula of [15]. In order to prove this we write the system as (9), and we have where the operator is a sectorial operator in with domain and has compact resolvent. In this context, the operator must be understood in the variational sense; that is, for every , and coincides with the fractional space of exponent [15]. Hereafter we denote by the norm on the space . Now, if we prove that the nonlinearity is well defined and is Lipschitz and bounded on bounded sets, we obtain the local existence for the initial data in .
Using and being locally Lipschitz together with and , we will prove the nonlinear terms, , and and satisfy , , and ; that is, is well defined, Lipschitz, and bounded on bounded sets. It is possible to prove this by considering . In order to prove these properties of the nonlinearity , let , and we note that where and adding ), and in , we have and adding in , we get and from the previous hypothesis on function , there exists such that and the rest is obvious.
Therefore, using the techniques of variations of constants formula of [15], we obtain the unique local solutions of (8) which are given by
with , and using again the results of [15], we get the regularity of solutions. In fact, from the smoothing effect of the equations, we have and , for some positive and any . Now, for we have , and since is well defined, Lipschitz, and bounded on bounded sets, we have and . Since is arbitrary, we obtain the regularity of the local solution.
Step 2. Now, we prove the solutions of (8) for every time .
To prove the global existence, we must show that the solutions are bounded in norm on finite time intervals. First, to obtain that the norm of is bounded in finite time, we note that multiplying the equations for the temperature by in and integrating by parts, we have: since .
Using Cauchy-Schwartz and the Young inequality and then the Poincaré inequality, since together with is the first nonzero eigenvalue of in , we obtain and using , we get and we conclude that the norm of in remains bounded in finite time.
Now, we note that differentiating the second equation of (5) with respect to , we obtain the same equations for , and considering now , we obtain
Thus, we show that the norm of in remains bounded in finite time. Then, using bounded for finite time, we prove that and remain bounded in finite time and we conclude.

2.1.2. The Case with No Diffusion:

The system now reads where ; that is, we consider Newton’s linear cooling law as in [68, 10], and it is no longer of a parabolic type system and is given by

To prove the system is well posed, we use the techniques from [16] considering the same transport equation for temperature in different thermosyphon models as [68, 10, 16].

We note that if is a given continuous function, then the equation for the temperature can be integrated along characteristics to obtain and plugging this into the nonlocal differential equation for the acceleration and into the equation for the velocity yields

We note that for and since in this space the translations are continuous isometries, (27) defines a continuous function of time with values in this space. Although we restrict ourselves to , many other choices of space are possible for solving problem (26). In fact any Banach space of 1-periodic functions of having zero mean and in which translations are continuous isometries can be used as an “admissible space” ; see [16]. In particular are admissible spaces between others. Then, we can prove Lemma 2.

Lemma 2. Let , fix , and assume that where is an “admissible space”; see [16], in particular . Then, the function given in (27), , is an integral solution of the PDE which is satisfied only if and are differentiable. In particular, if , then is continuous with values in and satisfies the PDE as an equality in , a.e. in time. Moreover, (27) satisfies the following properties:(i)(ii) if there exist positive constants and such that satisfy for all , then satisfies positive constant independent on time;(iii) we assume that and there exist positive constants and such that satisfy for all . If we also assume that are continuous in , then is a positive constant, and .

Proof. See [68, 10, 16]. Then, we have Theorem 3.

Theorem 3. Assume that is locally Lipschitz, , , and . Then there exists a unique solution of (26) satisfying and satisfies the PDE in the sense of (27). Moreover .

Proof. As noted earlier, we need to solve the fixed point problem on a space of continuous functions. More precisely, we take , endowed with the sup norm, with and to be chosen and prove that is a contraction on .
From (29) in Lemma 2 for , we have , and this, together with the local Lipschitz property of shows that for fixed , if is sufficiently small.
To show that is a contraction, it is clear that we must prove some Lipschitz dependence on with respect to .
First, we note that from , then verify (30), that is, for all with , and given , again from (31) in Lemma 2, we have With these we find that is Lipschitz on with a Lipschitz constant depending on and that tends to zero as and then is a contraction for small enough . Therefore, local well-posedness follows.
To prove the global existence, it is sufficient to prove that is bounded on finite time intervals, since from we find that is of Cauchy type as for finite . Consequently, the limit of exists in and the solution can be prolonged.
But again from (29) together with (18) and (19), we obtain boundedness on finite time intervals and global existence follows.
As noted earlier, if , then satisfies the PDE equation as an equality in a.e. in time. In particular, we have .

2.2. Boundedness of the Solutions and Global Attractor

In order to obtain asymptotic bounds on the solutions as , we consider the friction function satisfying the hypotheses from the previous section and we also assume that there exists a constant such that We make use of L’Hopital’s lemma proved in [11] to prove several results in this section.

Lemma 4 (L’Hopital’s lemma). Assume that and are real differentiable functions on on and .(i)If , then .(ii)If , then .

With this, we have Lemma 5.

Lemma 5. If one assumes that and satisfy the hypothesis from Theorem 1 or Theorem 3 together with (36), then with a positive constant such that if .

Proof. Integrating by parts we have and using Lemma 4 we obtain and from (36) we conclude.

Remark 6. We note that the conditions (36) are satisfied for all friction functions considered in the previous works, that is, the thermosyphon models, where is constant or linear or quadratic law. Moreover, the conditions (36) are also true for , as .

Now, we use the asymptotic bounded for temperature to obtain the asymptotic bounded for the velocity and the acceleration functions.

Theorem 7. Under the previous notations and hypothesis of Theorem 1 or Theorem 3, if one assumes also that satisfies (37) for some constant , then one has the following.
Part (I). General case:(i)In particular: If , then (ii) If and , then Part (II). If and also assume that there exists a positive constant such that , then for any solution of (5) in the space , one has(i)(ii)(iii) if , In particular, (5) has a global compact and connected attractor, , in .

Proof.
Part (I). General case.
(i) From (8) we have that and satisfies where . First, we rewrite (48) as with Next, for any there exists such that for any and integrating with , we obtain Using L’Hopital’s Lemma 4 proved in [11], we get Moreover, using again the L’Hopital’s Lemma 4 proved in [11], we get and from (51) together with (37), we conclude that for any .
(ii) From (47) together with Gronwall’s lemma, we get where . Consequently, for any , there exists such that for any , that is, for any , and using the previous results (i) we get (42).
Part (II). (i) From (23) together with (24), we get and by elementary integration we obtain (44). Using Part (I), the rest (ii) and (iii) are obvious. Since the sectorial operator defined in Section 2.1.1 has compact resolvent, the rest follows from [17, Theorems and ].

Remark 8. First, we note that the hypothesis about the function in Theorem 7, is satisfied when we consider Newton’s linear cooling law , where is a positive quantity; that is, , as [18]. Moreover, this condition is also satisfied if we consider where 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 (8), when we consider the general Newton’s linear cooling law without the additional previous hypothesis on ; but we assume that the friction function always satisfies (36).
In order to get this, we consider the Fourier expansions and observe the dynamics of each coefficient of Fourier expansions to improve the asymptotic bounded of temperature. In particular, we will prove for every locally Lipschitz and positive function and also for every (see (72) in Proposition 9) and for every friction function satisfying (36).

3. Asymptotic Behavior: Reduction to Finite-Dimensional Systems

We take a close look at the dynamics of (5) by considering the Fourier expansions of each function and observing the dynamics of each Fourier mode. Assume that and are given by the following Fourier expansions: while the initial data is given by .

Assume that is given by Then, we find that the coefficient in (60) is a solution of Since all the functions involved are real, we have , and . Therefore, (5) is equivalent to the infinite system of ODEs consisting of (61) coupled with

The system of (5) reflects two of the main features: (i) the coupling between the modes enter only through the velocity, while diffusion acts as a linear damping term, and (ii) it is important to note in this model that we have also the nonlinear term given by Newton’s linear cooling law.

We note that the system (5) is equivalent to the system (8) for acceleration, velocity and temperature. It is equivalent to the following infinite system of ODEs:

In what follows, we will exploit this explicit equation for the temperature modes to analyze the asymptotic behavior of the system and to obtain the explicit low-dimensional models.

3.1. Attractors and Inertial Manifolds

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.

A similar explicit construction was given by Bloch and Titi in [19] for a nonlinear beam equation where the nonlinearity occurs only through the appearance of the norm of the unknown. A related construction was given by Stuart in [20] for a nonlocal reaction-diffusion equation.

In order to get the inertial manifold for this system, we first improve the bounds on acceleration, velocity, and temperature of the previous section for all situations with and general locally Lipschitz function under the hypotheses of Lemma 5 for some , in particular satisfying (36).

We will prove in Proposition 9 that we have always an upper bounded for the temperature in independent of the velocity and the function , considered in Newton’s linear cooling law and also independent of the diffusion coefficient. That is,

Proposition 9. Under the previous notations, for every solution of the system (5), , and for every , one has(i)(ii)and   a positive constant such that ;(iii) In particular, one has a global compact and connected attractor where are the upper bounds for acceleration and velocity as given in (67) and (66) and .

Proof. From (61), we have with Thus, we obtain and we get . Using Theorem 7 together with , we get From this upper bounded for the velocity, we also have the upper bound for the continuous positive function , and using Part (II) from Theorem 7 we conclude.

We note from the previous result that we have always the upper bound for and from Theorem 7 for the velocity. Therefore, we can consider the upper bound for the continuous positive function ; we note that and we prove in Proposition 10 the bound of solutions to show the influence of diffusion coefficient .

Proposition 10. Under the previous notations, for every solution of the system (5), , and for every , one has Moreover, one has that and   positive constant such that ;

Proof. Using again if we assume that , then we obtain and we get that is, if , then .
In particular for every and also ; that is, .
Then, we also have that
Finally, we note that if , then given , there exists such that for every , and integrating in , we obtain , for any , and Using again Theorem 7, we conclude.

We note that if , then . If we observe in the numerical experiments, as is bigger, that the solution becomes stable or periodic.

As a consequence, we have the following result on the smoothness of the attractor of (5).

Corollary 11. (i) If , then for every .
(ii) If is the global attractor in the space , then for every , with , one gets In particular, if with , the global attractor and is compact in this space.

Proof. (i) From (77) we have , therefore, if , then for every and .
(ii) We note that from (i), if , then , and therefore ; that is, where are the upper bounds for acceleration and velocity as given in (74) and (73). But, the set is compact in since for any sequence in we can extract a subsequence that we still denote such that it converges weakly to a function and such that for any , the Fourier coefficients verify as , where is the th Fourier coefficient of . Therefore, and for every integer , where denotes the norm in . Hence, the first term goes to zero as and the second term can be made arbitrarily small as . Consequently, and in , and the result follows.

Note that this result reveals in particular the asymptotic smoothing of (5). In the next result we will prove that the dynamical system induced by (8) in the phase space , , has an inertial manifold. According to [21], we have the following definition.

Definition 12. Let , be a nonlinear semigroup in a Banach space in that has a global attractor . Then, a smooth manifold is called an inertial manifold if(i) is positively invariant, that is, for every ;(ii) contains the attractor, that is, ;(iii) is exponentially attracting in the sense that there exists a constant such that for every bounded set there exists such that for every .

See, for example, [21, 22].

Assume the ambient temperature given by where , that is, with for every with , since .

Then, we denote by the closed linear subspace of spanned by and consider the following spectral decomposition in , where denotes the projection of onto and the projection onto the space generated by . Note that (8) is equivalent to

Note that from (61), if , then the th mode for the temperature is damped out exponentially, and therefore the space attracts the dynamics for the temperature. This is precisely stated in the following result.

Theorem 13. Assume that and . Then, the set is an inertial manifold for the flow of in the space . Moreover, if the inertial manifold has the exponential tracking property; that is, for every , there exists such that if , are the corresponding solutions of (8), then in . In particular if is a finite set, the dimension of is , where is the number of elements in .

Proof. In order to prove that is invariant, we note if , then , and therefore , from (70), we get that , for every ; that is, . Thus, if , then for every , that is, invariant manifold.
We consider the decomposition in , , where is the projection of on and is the projection of on the subspace generated by ; that is, and .
From (77) taking into account that for , we have that for every implies that there exist positive constants , such that , that is, in if .
Therefore, we have in particular that as with exponential decay rate . Thus, also attracts with exponential rate , since
To prove the exponential tracking property just note that the flow inside is given by setting where ; that is, Therefore, if , then , and given and , the solution of (84), we decompose and . Then, we consider and it is still a solution of (85). Hence, and the right-hand side is of order . In particular, if the set is finite, then the inertial manifold is of finite dimension and the flow inside is equivalent to the finite system of ODEs given by (85). Thus, the theorem is proved.

3.2. The Reduced Subsystem

Under the hypotheses and notations of Theorem 13, we suppose that with for every . Then, . So, the evolution of velocity and acceleration depend only on the coefficients of which belong to the set . Note that in (85) the set of equations for with together with the equation for and is a subsystem of coupled equations denoted by the reduced subsystem.

Thus, we will reduce the asymptotic behavior of the initial system (5) to the dynamics of the reduced explicit system (87) when we consider the relevant modes of temperature . Consider where , , and as we consider only the real functions. After solving this, we must solve the equations for which are linear autonomous equations.

Now, we will show the modes in that will play an essential role in the dynamics. With the previous notations, we further decompose as follows: where is the projection onto the space generated by and is the projection onto the space generated by . We denote by the projection and . With these notations and decomposing as , (8) and (84) can be decomposed as a system of the form Since and setting , the first four equations give the flow inside the inertial manifold ; that is, they are equivalent to (84) while the first three are the only nonlinearity coupled equations. Therefore, once this subsystem is solved, the other unknowns are determined through linear nonhomogeneous equations.

To make this idea more precise in terms of semigroup and attractors, we proceed as in [14]. We denote by the semigroup generated by (8) on and by its restriction to the inertial manifold , that is, the semigroup generated by (90). Consider

We will find a reduced semigroup on the reduced space , denoted as , that, in a sense, determines the asymptotic behavior of and therefore that of . Note that and have the same attractor, while the dimension of the space might be much smaller than that of . The next result states, in particular, that the attractor of the full system can be reconstructed from the attractor of the reduced one.

Proposition 14. With the previous notation, one has the following conditions.(i)The system of equationsdefines a nonlinear semigroup, denoted as , on that can be identified with restricted to .(ii) If denotes the maximal attractor of (8), then is the maximal attractor of (91). Moreover, where is continuous.(iii) If the set is finite, (91) is equivalent to a system of complex ODEs of the form (87). Consequently, the asymptotic behavior of (8) is described by an explicit system of ODEs in with an even number. In particular, if , and for every , one has that the associated solution verifies and in , where is the unique solution in of the equationMoreover, if , one gets and .

Proof. (i) Working as aforementioned, we prove that the semigroups and are well defined and prove the existence of attractor . Using the techniques of [16] we can work as in [14] to prove (ii). Then (iii) we note that , and since and , then the set is a symmetric set and has an even number of elements that we denote by . Therefore, the number of the positive elements of , , is .
Note that . Thus, the dynamics of the system depends only on the coefficients in . Moreover, the equations for are conjugates of the equations for , and therefore we have that Taking real and imaginary parts of , that is, employing real variables, , we have a system in with .
If , then from the equation for the velocity we have
Moreover, from the equation for the temperature in (5), we have that the function satisfies the equation We can multiply by in , and taking into account that , since is periodic, we have and using Cauchy-Schwartz and Young inequality with and then Poincaré inequality, since , together with , we have that Next, using , we prove that in .
Now, we multiply (97) by in . Integrating by parts, applying Young inequality and taking into account again that since is periodic, together with , we obtain that for every with . Thus, working as aforementiond and taking into account that , we get in ; that is,.

Next, we pay attention to the other modes for the temperature where . Also, note that these modes are determined as solution of the linear nonhomogeneous equations with initial data . Therefore, we call these the slave modes.

We will show in Proposition 15 that the dynamics of these modes are completely determined by the solution of (87), in the sense that the solution will have only one asymptotic behavior as time goes to infinity.

Proposition 15. Assume that is a solution of (87). Then, for any , there exists a solution of (101), denoted as , such that for every and for any other solution of (101), at an exponential rate independent of , as . Moreover, if , that is, if , then ; that is, this subset of the slave mode is damped out exponentially. In particular, if is a stationary or periodic (resp., quasiperiodic, almost periodic) solution, then can be chosen such that it is stationary or periodic with the same period (resp., quasiperiodic, almost periodic with a set of frequencies contained in those of .

Proof. Define as a solution of (101) with an initial condition satisfying . In particular, if , that is, if , then . Then for any other solution of (101), satisfies the homogeneous equation and which proves the statement.
If the solution of (87) is stationary, that is, independent of time, then we choose to be the solution of and the result follows.
If is a periodic solution of (87), then since the are the solutions of linear scalar differential equations of the form , with and period of the same period and the homogeneous equation is stable, the result follows from Fredholm’s alternative. For the quasiperiodic or almost periodic case, the result follows from Theorem  6.6 in [23].

Remark 16. Taking real and imaginary parts of , and as the asymptotic behavior of the system (5) is given by a reduced explicit system in , where , given by where , , and .
Observe that from the previous analysis, it is possible to design the geometry of circuit and/or the external heating by properly choosing the functions and/or the heat flux and the ambient temperature so that the resulting system has an arbitrary number of equations of the form .
Note that the set can be much smaller than the set , and therefore the reduced subsystem may possess far fewer degrees of freedom than the system on the inertial manifold. Also note that it may be the case that and are infinite sets, but their intersection is finite. Also, for a circular circuit, we have , that is, and then is either , or the empty set. Also, if in the original variable for (5) is constant, we get for any choice of .

The physical and mathematical implications of the resulting system of ODEs which describe the dynamics at the inertial manifold need to be analyzed numerically. The role of the parameter which contains the viscoelastic information of the fluid deserves special attention and will be the aim of the next section.

4. Numerical Experiments

In this section, we describe the results of the numerical experiments obtained using the MATHEMATICA package [24] for the resolution of the differential equations, using a fourth-order explicit Runge-Kutta method for stiffness equations following the method used in previous works [5, 15]. We will solve a system of ordinary differential equations which are the projection of the partial differential equations (5) on the inertial manifold derived in the preceding sections. All the variables and equations that we deal with are nondimensional. 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).

Specifically, we are integrating the system of (87) where we consider only the coefficients of temperature with (relevant modes). Then, where , and since all the physical observable are real functions. In particular, we will consider a thermosyphon with a circular geometry; so and . Consequently, we take and omit the equation for . Hence, where the unknowns are (the acceleration of the fluid), (velocity of the fluid), and (the Fourier mode of the temperature). More complex geometries will result in higher dimensional dynamics on the inertial manifold. It is beyond the scope of the present work.

In order to reduce the number of parameters we make the change of variables and we define the real and imaginary parts of the equations in the following way: with . Therefore, our central results correspond to the system of equations Note that it is a system of four equations with four unknowns where we need to make explicit choices for the constitutive laws for both the fluid-mechanical and thermal properties. Thus, for the friction law and heat flux we will take the ones used in [5, 15]. For the numerical experiments, which are of a similar model of thermosyphon for a fluid with one component they use the functions and . The function has a clear physical meaning; it interpolates between a low Reynolds number friction law (in which the overall friction is linear (Stokes friction law) and high Reynolds number (in which the friction is a quadratic law).

and refer in this model to the position-dependant () heat flux inside the loop and will be used as tuning parameters. Without loss of generality, we will assume that in order to simplify in analogy with Lorenz’s model, as it is shown in [5, 15] (changing and simultaneously only results in a change in the phase of initial temperature profile).

We have carried out two different sets of numerical experiments with regard to heat diffusion. The first numerical experiments are carried out keeping the heat diffusion to zero as it was done in [11]. And the second numerical experiments are performed with heat diffusion. The initial conditions are fixed as . This split would appear naive as diffusion tends to smooth the solution however, as the order of the equations changes in the presence of diffusion (from first to second order due to the Laplacian), it is worth studying both cases separately.

Numerical analysis has been carried out keeping the viscoelastic coefficient as the tuning parameter ranging from 100 to 0.0001 and the heat flux also as a tuning parameter ranging from 1 to 10000. The impact of on the system has been keenly observed for various parameters of time , as short as 50 time units and as long as 5000 time units. 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, that is, transient irregular behavior followed by convergence to equilibria, periodic behaviors, and chaotic progressions.

4.1. Dynamics of the Thermosyphon without Heat Diffusion ()

In this section, we summarize some of the outcomes of the model equations. As we have mentioned earlier, this behavior is highly sensitive to the choice of parameters. Thus, we present those results in different subsections accounting for the most relevant signature for each set of numerical experiments.

4.1.1. Chaotic Behavior of the Model for Large Values of

The simulations of the numerical experiments done for large values of , for instance, ranging from 2 to 1000, show that the system exhibits a pattern of chaotic behavior. For relatively large values of viscoelastic components, chaotic behaviors are observed. For all the values of heat flux , starting from 1 to 10000, this chaotic behavior is observed (see Table 2). To illustrate the previously stated chaotic behavior, we observe the plots obtained for the values of and .

In Figure 1, we show a time graph of acceleration for a large value of the viscoelastic parameter . The acceleration ranges from −3.5 to 3.5. Since the very beginning, it displays a chaotic behavior. The curve is not very erratic although it does not show any sort of periodicity. As velocity is the time integral of acceleration, in Figure 2, the curve does not present abrupt changes close to some maxima and minima, but the nonperiodic features are also captured by this observable.

In Figure 3, we show a phase-diagram plot for the real and imaginary parts of the temperature. As expected, it also exhibits a nonperiodic pattern in which the trajectory in this phase-plane moves inwards and outwards the graph. This graph illustrates the complex underlying dynamics of the attractor (of which Figure 3 is a two-dimensional projection).

This sort of behavior remains similar for other values of as long as the viscoelastic parameter takes larger than unity values as summarized in Table 2. In other words, the elastic effects introduce a memory effect in the dynamics which avoids the system to fully stabilize but, rather, viscoelasticity sustains the chaotic pattern. This memory effect can be understood from (3).

To sum up this section, large values of the viscoelastic parameter result on sustained chaotic behaviors. The dynamics becomes more complex, characterized in all the cases by periods of chaos and of violent oscillations, giving an idea of the complexity of the solutions of the system under these variables. In the detailed analysis of the evolution of the acceleration, velocity, and temperature, we say that the chaotic behavior of the system reveals the chaotic nature of the viscoelastic fluids.

4.1.2. Transient Irregular Behavior Followed by Stable Behavior for to 1

For values of ranging from to , the system tends towards a stable fixed point, although it is still chaotic in the initial stages. This transient irregular behavior followed by equilibrium is shown in Figures 4, 5, and 6.

The general behavior of acceleration is that it has a chaotic outburst in the initial stages. But as time progresses it tends to stabilize, attaining equilibria. The velocity too, in the initial stages, when the time period is less than 20 units, is very inconsistent and at times unpredictable. But as the time progresses, velocity converges to a stable fixed point. Interestingly, this fixed point for the velocity is not trivial () but, on the contrary, it depends strongly on the choice of the parameters. Specifically, in Table 1, we summarize this asymptotic value. Physically, this means that the fluid inside the thermosyphon performs sustained flow in the same direction over time. This stage could be identified by the existence of convective rolls in a fully spatially extended system.

It is worth noting in Table 1 that both columns have the same value for the asymptotic velocity. This is a signature that viscoelastic effects do not play any role in this case (so memory effects are damped out after the early chaotic transient).

Another comment regarding Table 1 concerns the role of the parameter . Roughly, accounts for the scale of temperature gradients inside the thermosyphon. These results suggest that higher temperature gradients produce higher values of the sustained stationary velocity. Actually, as shown in Figure 6, the equilibrium velocity scales with nontrivially (as a power law with exponent approximately).

To sum up this section, we conclude that although the system has a chaotic initial transient, it tends to stabilize at longer times reaching a temperature gradient dependent equilibrium velocity. Notwithstanding, this asymptotic velocity depends nontrivially on temperature as a power law, being this a signature of the underlying nonlinearity of the equations.

4.1.3. Transition to Periodic Pattern of Behaviors for Small Values of

When the viscoelastic effects are gradually less important (values of between 0.01 and 0.0001), the system exhibits different behaviors as a function of the temperature gradient (see Table 2). For lower values of , the system behaves in a similar fashion as for higher values of . On the contrary, for larger values of , the system displays a periodic pattern. As in the previous case, both long-term behaviors may be preceded by an initial chaotic transient (probably caused by viscoelasticity).

To illustrate this, in Figure 7, we show a periodic (nontrivial) behavior for . In some cases, the initial transient cannot be distinguished from the periodic one and we refer it simply to periodic type in Table 2.

In Figure 8, we show the phase-space of the complex components of the temperature. Although at first sight it resembles a typical chaotic attractor motif, after the mentioned initial transient, the trajectories are overlapped for longer times.

Thus, in order to summarize the information covered in the last three subsections, we collect all the outcomes of the model in Table 2.

4.2. Dynamics of the Thermosyphon with Heat Diffusion ()

In this case, to avoid unnecessary repetitions in the text, we focus on the main differences between this case and that in Section 4.1.

Thus, in this second set of numerical experiments we introduce a nonzero value for the thermal diffusivity, . The role of diffusion is to reduce temperature gradients. This can be seen in a hand-waving way by realizing that namely, the Laplazian can be understood as the flux of temperature created by a temperature current, . This current is larger in those regions where the temperature variations are also larger. Thus, the system tends to reduce those differences. As shown in Table 3, the variety of the behaviors is clearly less rich than in the case when .

So, here we will only illustrate the most interesting behaviors with two examples. The first one takes the values of heat diffusion , , and . As shown in Figure 9, the acceleration performs a series of damped oscillations that eventually stabilizes.

Similarly, the second example (Figure 10) takes the values , , and . The behavior is qualitatively equal but the period of the oscillations is enlarged. In Table 3, we summarize the interaction between the tuning parameters.

To sum up this section, we have found that greater values of the heat diffusion, , smoothen the dynamics of the system which, invariably, tends to stabilize, either reaching an equilibrium steady state or stable periodic orbits. These two behaviors are governed by the value of the temperature gradient which in a similar fashion as in Section 4.1 tends to produce richer behaviors for large values.

5. Conclusion

In this work, we have derived a novel system of equations to study the behavior of a viscoelastic material inside a thermosyphon. This model serves as a preliminary simplification of a more complex fully spatially extend system. This model illustrates the presence/absence of complex chaotic behaviors and also enables to relate them with the underlying viscoelastic (memory effects).

The main result is that we are able to prove that the original system (which involves both ordinary and partial different equations) possesses an inertial manifold in which the dynamics can be accurately described by a system of ODEs. By numerical integration of the reduced equations, we have been able to better understand the role of viscoelasticity (as opposed to a simpler Newtonian fluid) through the parameter . This parameter is a nondimensional version of the so-called Maxwellian viscoelastic time [4] which accounts for memory effects.

Our results suggest that in the absence of heat diffusion (), a large value of drives the dynamics to chaotic behaviors for all the physical observable (acceleration, velocity, and temperature). As the value of gradually decreases, the system is no longer chaotic but stable or periodic. Notably, these results cannot be understood in terms of a boundary layer theory as the attractor of the dynamics changes dramatically when the second-order derivative term is introduced in (5).

Physically, this induction of chaotic behaviors can be rationalized in terms of the memory effects inherent to viscoelastic models explicitly shown in (3). Thus, in the same way as delayed equations are known to produce chaos, even in the simplest situations, viscoelasticity apparently produces the same kind of transition (see, e.g., [25]).

Other interesting results are related to the effect of heat diffusion. We have found that as the heat diffusion is not zero, the system tends to stabilize either to a fixed equilibrium point or to a (-dependent periodicity) periodic orbit. This abrupt change in the qualitative behavior of the system is, again, due to the fact that the term is a singular perturbation that changes the type of equation for the temperature from hyperbolic to parabolic, which justifies the explicit analysis of each situation separately in Section 2.

Our work can be generalized in many different ways, from changing the constitutive equation (from Maxwellian to other more complex situations) or to include shear-thinning effects [4] common to many non-Newtonian materials. Shear thinning is the manifestation of a shear-rate-dependent viscosity. Thus, it is commonly observed that many fluids reduce their resistance to flow for large enough imposed stresses (in our case, temperature gradients), for instance, tooth paste, paint, or lava. In principle, shear thinning means that the viscosity is smaller for higher velocities what could help to sustain oscillations once they have been formed inside the loop. However, as shown in the present work, this sort of hand waving analysis might lead to erroneous conclusions, and a detailed consideration of those effects deserves further (and deeper) investigation. These generalizations will be the aim of future works.

Appendix

Boundary Layer Theory

The system of (5) can be seen as a singular perturbation problem provided that when , the order of the differential equation reduces to one. In order to provide some insight about the solutions of those equations, a standard boundary layer theory is customary (see, e.g., [26]).

In particular, one splits the problem into two parts: the inner problem and the outer problem. The inner problem is defined as the dynamics of the system for times up to in which the term is dominant. This regime is strongly dominated by transient impulsive changes in the physical variables. Following [26], we define a new time scale . So, the first equation in (5), up to order , is given by whose solution is with a constant that can be determined by matching with the outer solution. For instance, if one assumes, the velocity has initial condition and the matching time is , then it is straightforward to see that the velocity at a matching time converges exponentially to .

Besides, the outer problem is defined as the naive approximation . Thus, the system of (5) reduces to that in [11] with an effective initial condition given by . So one would expect the same qualitative behaviors as in that work.

For instance, if the parameter , the model in [11] predicts a stable behavior. However, as shown in Table 2 (second row), increasing the value of induces a chaotic behavior. Naturally, different values of the parameters (including ) would provide different values of for which the trivial case cannot account for the observed dynamics (even for small enough values of ).

This simple analysis shows the intrinsic complexity of the physical problem when viscoelasticity is considered and, more importantly, the need to study every parameter set in detail to provide an accurate description of the type of dynamics in which the system evolves.

Acknowledgments

This work is partially supported by projects MTM2009-07540, GR58/08 Grupo 920894 BSCH-UCM, Grupo de Investigación CADEDIF, and FIS2009-12964-C05-03, Spain, and partially supported by Grant MTM2012-31298 from Ministerio de Economia y Competitividad, Spain.