An analysis of the relaxation oscillations of a nonlinear thermosyphon

Author’s Accepted Manuscript An analysis of the relaxation oscillations of a nonlinear thermosyphon Moundheur Zarroug, Peter Lundberg, Fariba Bahrami www.elsevier.com/locate/nlm

PII: DOI: Reference:

S0020-7462(16)30133-0 http://dx.doi.org/10.1016/j.ijnonlinmec.2016.09.003 NLM2703

To appear in: International Journal of Non-Linear Mechanics Revised date: 2 September 2016 Accepted date: 8 Cite this article as: Moundheur Zarroug, Peter Lundberg and Fariba Bahrami, An analysis of the relaxation oscillations of a nonlinear thermosyphon, International Journal of Non-Linear Mechanics, http://dx.doi.org/10.1016/j.ijnonlinmec.2016.09.003 This is a PDF file of an unedited manuscript that has been accepted for publication. As a service to our customers we are providing this early version of the manuscript. The manuscript will undergo copyediting, typesetting, and review of the resulting galley proof before it is published in its final citable form. Please note that during the production process errors may be discovered which could affect the content, and all legal disclaimers that apply to the journal pertain.

An analysis of the relaxation oscillations of a nonlinear thermosyphon Moundheur Zarrouga,∗, Peter Lundberga , Fariba Bahramib a

Department of Meteorology/ Physical Oceanography, Stockholm University, 106 91 Stockholm, Sweden b Department of Applied Mathematics, Faculty of Mathematics, University of Tabriz, Tabriz, I. R. of Iran

Abstract The oscillatory behavior of an asymmetrically forced thermosyphon constituted by two connected vessels has been subjected to an asymptotically valid analysis using the vessel-volume ratio as expansion parameter. Due to the structure of the governing equations, the problem could not be dealt with using standard techniques; instead a phase-plane analysis was conducted. The analytically determined corrections to the previously established lowest-order discontinuous results proved to be useful even for comparatively large values of the expansion parameter. The relationship between these asymptotically valid corrections and the physics underlying the relaxation oscillation as well as the behavior of the system for strong thermal forcing is discussed. The study is concluded by an overview of some specific inconsistencies associated with the discontinuous lowest-order analysis and how these were alleviated by the asymptotically valid corrections. Keywords: Thermosyphon, Nonlinear dynamics, Relaxation oscillations, Asymptotic analysis 1. Introduction A thermosyphon is basically a heat-transfer device employing either standard thermal convection (single-phase devices) or the principle of evaporation and condensation of the working fluid (two phases). A single-phase ∗

Principal corresponding author Email address: [email protected] (Moundheur Zarroug)

Preprint submitted to International Journal of Non-Linear Mechanics September 9, 2016

thermosyphon is constituted by a closed loop containing the working fluid which is subjected to heating and cooling in its lower and upper regions, respectively. Practical realisations of two-phase thermosyphons commonly are based on a rod-like vessel partially filled with the working fluid. When this device is heated at the bottom and cooled at the top, the cycle is closed by the return flow of the condensate taking place along the interior walls of the closed tube. Engineering applications of thermosyphons began more than a century ago when this principle was utilized for cooling systems aimed at maintaining internal-combustion engines at optimal working temperatures. In these early days empirical experience formed the basis underlying system design. Systematic theoretical and experimental investigations of thermosyphons were initiated in the 1950s when it was realized that these had an essential role to play when designing the cooling systems of light-water as well liquidmetal breeder nuclear reactors [1]. Research dealing with these matters has been carried out until the present day, see e.g. the review in [2] and the experimental results reported in [3], [4]. In the wake of the oil embargo in the early 1970s, a renewed interest in solar-heating systems arose. Since these devices almost invariably are based on single-phase closed-loop devices with water as the working fluid, this led to a considerable body of theoretical as well as experimental work on thermosyphons being undertaken, cf. e.g. [5] and [6], [7], respectively. More recently the exploitation of natural resources in permafrost regions in conjunction with expectations of climate change has led to the widespread use of rod-like two-phase thermosyphons being used to ensure stability of the foundations (viz. by preventing permafrost degradation) when constructing pipelines, roads and railways, telecommunication towers, dwellings etc. in Arctic and Subarctic environments. A useful review of these matters is provided by [7], and in [8] a recent laboratory investigation with practical implications is reported. The cooling of electronic equipment (power electronics as well as the integrated circuits and processors used in computer applications) has, not least with the advent of miniaturization, assumed increased importance, and the use of two-phase thermosyphons to maintain adequate working temperatures of the devices in question has become widespread. Since the theoretical background of thermosyphon operation was regarded as comparatively wellknown, the main research efforts in this field have been devoted towards experimental investigations, see e.g. [9], [10]. 2

Although the brief review above has focused on the engineering use of thermosyphons, it should also be underlined that a large body of work on closed-loop single-phase thermosyphons has been undertaken from a purely fluid-dynamical perspective. This has partly been done in the hope that these studies would be capable of clarifying the transitions from ordinary B´enard convection to fully developed turbulence as the thermal forcing gradually is increased, but also to examine double-diffusive convection. The first investigation of this type was conducted by Welander [12], and has been followed by many others (see e.g. [13]-[15]), frequently invoking insights concerning the route towards chaos based on modern dynamical-systems theory. From the references provided above it may be recognized that thermosyphons have been investigated using a variety of methods, namely analytically, numerically and experimentally. Despite the apparent simplicity of these devices a number of difficulties remain to be investigated, subtleties stemming from the fundamentally nonlinear nature of the equations describing these thermal systems. The present analytical study will deal with the oscillatory behaviour of a closed-loop thermosyphon based on the approximately quadratic equation of state for fresh water at temperatures around 40 C (viz. the environment commonly encountered in northerly regions by heat pumps feeding on lakebottom water). This working fluid is contained in two equally high, wellmixed truncated cylinders of volumes V1 and V2 , separated by a thermally conducting common wall with connecting tubes at top and bottom, cf. Fig. 1. Differential thermal forcing via the external walls of the vessels yields a density difference which drives a flow through the lower tube, the upper acting as a passive conduit. Under certain parametric conditions (viz. slightly asymmetric forcing around the temperature of maximum density and a large volume difference between the vessels) the resulting oscillatory behavior of the fluid temperatures T1 and T2 in the vessels of volumes V1 and V2 , respectively, is of a relaxation character, cf. [16]. The time evolution of T1 and T2 (representing deviations from 40 C) is governed by the following non-dimensionalized ordinary differential equations: δ

dT1 θ = Ra | T22 − T12 | (T2 − T1 ) + ( − T1 ) + κ(T2 − T1 ) = F (T1 , T2 ), (1) dt 1−θ

dT2 1 = −Ra | T22 − T12 | (T2 − T1 ) + ( − T2 ) − κ(T2 − T1 ) = G(T1 , T2 ). (2) dt 1−θ 3

Note that here Ra represents a slight modification of its counterpart as previously defined by the present authors [17], since a redundant parameter has been absorbed by not only the temporal variable t, but also by the Rayleightype number Ra; these changes do not involve any alterations of the overall characteristics of the physical system. In both of these governing equations the first r.h.s. term represents the advective heat flux to/from the vessel, the second the externally forced conductive heat flux affecting the vessel in question, and the third the conductive heat flux through the common separating wall of the system. In Eq. (1) the ratio between the vessel volumes δ = V1 /V2 , is taken to be small, cf. the relaxation character of the problem. Ra is related to the strength of the advective heat exchange between the vessels through the connecting tubes. θ =T1ex /T2ex is the ratio between the external forcing temperatures giving rise to the Newtonian heat fluxes affecting the working fluid through the external walls of the two vessels, and κ is the quotient between the internal and external coefficients of thermal conductivity characterizing the shared boundary and the outer walls of the vessels, respectively. The present investigation solely deals with periodic behavior of the system, and θ as well as κ will be kept fixed and essentially play passive roles in the analysis. (Their values must, however, be chosen such that an oscillatory state is realized; note in particular that for δ < 1, the external thermal forcing must be such that −1 < θ < 0.) The study will, however, take into account the effects of varying the Rayleigh-type number Ra and the volume ratio δ. It is important to underline that this oscillator is fundamentally nonlinear in the sense that if the working fluid is taken to have a linear equation of state (i.e. were its density to be monotonous function of temperature), the thermosyphon would not display any oscillatory behavior but rather manifest a stationary circulation. A consequence is that Eqs. (1) and (2) after linearization of this type do not reduce to a linear oscillator which can be used as the basis for a perturbative analysis, this in marked contrast to the van der Pol [16] equation. A further problem is that the right-hand members of Eqs. (1) and (2) are non-analytical, which makes it necessary to deal with the problem piecewise in the two sectors T22 > T12 and T22 < T12 . This precludes reformulating the problem in terms of a second-order ODE in one of the temperature variables. When attempting to pursue an asymptotic analysis of the problem it is hence not feasible to apply the established techniques recently discussed by e.g. Grasman [18] and Verhulst [19]. In what follows the problem will instead be 4

analyzed in the (T1 , T2 ) phase-plane, cf. Mishchenko and Rozov [20]. Along these lines two lowest-order analyses of the thermosyphon have previously been carried out (cf. Lundberg et al. [17], Bahrami et al. [21]), both of which imposed physical constraints on the system. In [17] it was assumed that δ << 1, corresponding to the larger-vessel temperature remaining constant during the rapid phases of the oscillation. Some pertinent results illustrating this investigation are summarized in Fig. 2, showing that in the (T1 , T2 )-plane the relaxation oscillation to lowest order can be seen as constituted by periods of slow motion interspersed by rapid jumps (discontinuities) from one branch of F (T1 , T2 ) = 0 to the other. Smaller values of δ correspond to shorter time-scales for the jumps, which in the limit of δ → 0 take place parallel to the T1 -axis. In its motion along the two branches of F (T1 , T2 ) = 0 the phase-point loses stability at the cusp A, where these branches merge, and at the junction point C, corresponding to the minimum of T2 with respect to T1 on F (T22 > T12 ) = 0. In both cases rapid motion ensues, whereafter the phase-point resumes its stable motion on the other branch of F (T1 , T2 ) = 0 at the drop points B and D, respectively. In [21] another physical constraint was imposed, namely assuming that the total heat content of the system, T1 + δT2 , remained unchanged during the rapid phases of the oscillation, i.e. the jumps were taken to be adiabatic. In [17] and [21] additional information of physical nature was introduced to handle the problem of discontinuities. The aim of the present study is to relax these assumptions and instead deal with the relaxation-oscillation problem from a purely mathematical standpoint by applying asymptotic analysis to Eqs. (1) and (2). This will be accomplished by deriving asymptotically valid corrections of order delta to the lowest-order results reported in [17] and illustrated by Fig. 2. As shown in Appendix A, it is a straightforward matter to resolve the algebraic lowest-order equations yielding the precise locations of the cusp A and the junction point C as well as the drop points B and D. Once this has been accomplished, the derivation of the asymptotically valid corrections commences by examining the solution characteristics during the slow motion in the phase-plane sector T22 > T12 . Solution behavior in the region of the lowest-order junction point C is examined in Section 3, whereafter the subsequent rapid motion and behavior adjacent to the drop point D are dealt with in Sections 4 and 5. The slow motion in the sector T22 < T12 terminated by a ”jump” towards the analogous segment in the sector T22 > T12 is analyzed in Sections 6 and 7, and solution behavior in the region of the 5

drop point B is examined in Section 8. This completes the formal analysis of the first-order corrections to be applied to the lowest-order approximate solution. In Section 9 the resulting approximated limit cycles are compared to the lowest-order as well as numerical results [17]. The physics underlying the various stages of the relaxation oscillation are discussed in Section 10, and Section 11 shows how and why the asymptotically valid representations of the limit cycle are improved for stronger forcing, i.e. larger values of Ra. Some alternative considerations of a thermodynamic nature are applied to the oscillation in Section 12. The study is concluded by Section 13 focusing on the shortcomings of the lowest-order discontinuous solution reported in [17] and [21] and a discussion of how these were alleviated by asymptotically valid corrections. 2. Slow motion in the sector T22 > T12 To deal with solution behavior near the two branches of F (T1 , T2 ) = 0, a unified form of the governing equations in this phase-plane sector is obtained by dividing the pertinent version of Eq. (2) with that of Eq. (1): 1 dT2 ∓Ra(T22 − T12 )(T2 − T1 ) + 1/(1 − θ) − T2 − κ(T2 − T1 ) , (3) = Hν (T1 , T2 ) = δ dT1 ±Ra(T22 − T12 )(T2 − T1 ) + θ/(1 − θ) − T1 + κ(T2 − T1 ) where the upper signs corresponding to ν = 1 pertain to the phase-plane sector T22 > T12 , the lower signs and ν = 2 to T12 > T22 . By introducing the power series (0)

(1)

T2 (T1 , δ) = T2 (T1 ) + δT2 (T1 ) + ... ,

(4)

Eq. (3) can be reformulated as (0)

(1)

dT KT + δLT + ... dT2 . + δ 2 + ... = δ dT1 dT1 KN + δLN + δ 2 MN + ... In the sector T22 > T12 of present interest, 1 (0) 2 (0) (0) (0) − T2 ) − κ(T2 − T1 ), KT = −Ra(T2 − T12 )(T2 − T1 ) + ( 1 − θ   (1) (0) 2 (0) LT = −T2 Ra(3T2 − T12 − 2T1 T2 ) + (κ + 1) , (0) 2

KN = Ra(T2

(0)

(0)

− T12 )(T2 − T1 ) + θ/(1 − θ) − T1 + κ(T2 − T1 ), 6

(5)



 (0) − T12 − 2T1 T2 ) + κ ,     (2) (0) 2 (0) (1) 2 (0) = T2 Ra(3T2 − T12 − 2T1 T2 ) + κ + T2 . Ra(3T2 − T2 ) (6) (1)

LN = T 2 MN

(0) 2

Ra(3T2

Reordering Eq. (5) in powers of δ yields: δ0 :

(0)

KN dT

dT2 = 0, dT1

δ1 :

(0)

LN

dT2 = KT , dT1

(1)

δ2 :

LN

(0)

dT2 dT + MN 2 = LT . (7) dT1 dT1

(0)

Since dT21 = 0 during the slow phases, the lowest-order order result reduces to KN = 0. This relationship corresponds to the slow-motion branch F (T22 > T12 ) = 0 used for the discontinuous analysis, cf. the previous section. By simplifying KT using the lowest-order result KN = 0 as well as this equation differentiated with respect to T1 , viz. (0) 2

(0)

(0)

dT2 Ra(T2 + 2T1 T2 − 3T12 ) + κ + 1 , = (0) 2 (0) dT1 Ra(3T2 − T12 − 2T1 T2 ) + κ

(8)

the first-order correction term was found to be (1) T2 (T1 )

(0)

(0)

(1 + θ)/(1 − θ) − T1 − T2

=

(0) 2

(0)

Ra(2T1 T2 + T2

(0)C

− 3T12 ) + κ + 1

.

(9)

(0)C

As T2 → T2 and T1 → T1 the numerator remains non-zero whereas the denominator approaches zero, this in accordance with the relationship (0)C (0)C derived in Appendix A between T1 and T2 : (0)C T1

(0)C

T = 2 3

1 − 3



(0)C 2 4(T2 )

κ+1 +3 Ra

1/2 ,

and consequently the correction term becomes unbounded in a neighborhood of C, and another expansion is required near this junction point. Observe that for notational convenience, the coordinates of the lowestorder junction and drop points A, B, C, and D will henceforth (comprising the appendices) not include the superscript (0).

7

3. Analysis adjacent to the junction point C The general thrust of the investigation of solution behavior in this junctionpoint region adheres to the approaches pioneered by Haag [22] and Dorodnicyn [23] for dealing with the van der Pol equation. The analysis near C, i.e (T1C , T2C ), hinges on this lowest-order junction point being a local minimum of T2 with respect to T1 on F (T22 > T12 ) = 0 and hence is characterized by dT1 /dT2 |T1C ,T2C = 0. Since ∂F/∂T2 = 2RaT2 (T2 − T1 ) + Ra(T22 − T12 ) + κ is positive definite in the neighborhood of C, the implicit-function theorem is applicable and yields that here dT1 ∂F ∂F =− / , dT2 ∂T1 ∂T2

(10)

which in turn implies that ∂F/∂T1 = 0. In this neighborhood F (T1 , T2 ) = 0 can be reformulated as T2 = τ2 (T1 ), where τ2 (T1C ) = T2C and F (T1 , τ2 (T1 )) = 0. Differentiating the latter identity twice yields 



τ2 (T1C ) = 0;



τ2 (T1C ) = −

FT 2 

1

F T2

= 0.

(11)

Consequently a Taylor expansion applied to T2 = τ2 (T1 ) assumes the form 

T2 =

T2C

1 FT12 C 2 − .  (T1 − T1 ) + ... 2 F T2

(12)

This relationship implies that in a sufficiently small neighborhood of C, the introduction of the new variables ξ and η:

ξ = where    2∂F/∂T2    M=  2 ∂ F/∂T 2 

1 T C ,T C 1 2

T1C − T1 ; M

η = T2C − T2 ,

    Ra(3T C 2 − 2T C T C − T C 2 ) + κ  2 1 2 1 =  ,   Ra(3T1C − T2C )

(13)

(14)

makes it feasible to here reformulate F (T22 > T12 ) = 0 as η = −ξ 2 , b1 ≤ ξ ≤ b2 , the boundaries of this interval to be determined when evaluating the limit 8

cycle. (Note that the reversed orientation of the coordinate axes serves a useful purpose since the forthcoming local analysis will be conducted in terms of an increasing ξ.) The resulting set of governing equations turns out to be δ

dξ 1 = − F (T1 (ξ), T2 (η)), dt M

(15)

dη = −G((T1 (ξ), T2 (η)), dt

(16)

and thus Eq. (3) can be recast as dη G (T1 (ξ), T2 (η)) =δ M. dξ F (T1 (ξ), T2 (η))

(17)

   F T1C − M ξ, T2C − η ≈ (η + ξ 2 )FT2 T1 (ξ), T2 (ξ 2 ) ,

(18)

Here

as recognized by the artifice of subjecting a slightly manipulated version of F (T1 (ξ), T2 (η)) to a Taylor expansion, viz.   F T1C − M ξ, T2C − η + ξ 2 − ξ 2 ≈ F T1C − M ξ, T2C + ξ 2 −   −(η + ξ 2 )FT2 T1 (ξ), T2 (ξ 2 ) + ... . (19) Here the first r.h.s. term is recognized as being equal to zero since η = −ξ 2 for b1 < ξ < b2 . The governing equations can thus be reformulated as dη γ(ξ, η) =δ 2 , dξ ξ +η

(20)

where γ(ξ, η) =

−M G(T1C − M ξ, T2C − η) . FT 2 (T1C − M ξ, T2C + ξ 2 )

(21)

Note that γ(0, 0) = γ > 0. Once this overall formalism for dealing with the region adjacent to the junction point C has been established, a detailed analysis of solution behavior here can be undertaken. 9

3.1. Analysis of the initial part of the junction section To deal with Eq. (20) the following expansion is introduced: η(ξ) = η0 (ξ) + δη1 (ξ) + ... .

(22)

Insertion of this expression and ordering in powers of δ yields: dη0 2 (ξ + η0 ) = 0, dξ dη0 dη1 2 (ξ + η0 ) = γ(ξ, η0 ) − η1 . : dξ dξ

δ0 :

(23)

δ1

(24)

In accordance with previous considerations it is now assumed that η0 = −ξ 2 , which yields η1 (ξ) = −

γ(ξ, −ξ 2 ) . 2ξ

(25)

This correction term becomes unbounded as ξ → 0, and hence cannot be used for the asymptotic analysis in the immediate neighborhood of the junction point C. Its validity is limited to the interval b1 < ξ < σ1 , where σ1 → 0 as δ → 0. In (T1 , T2 )-space the expansion of η(ξ) over the initial part of the junction section thus assumes the form   C T1 −T1 (T1C −T1 )2 , − M γ C 2 M M2 (T1 − T1 ) C T2 (T1 , δ) = T2 + +δ + ... .(26) M2 2(T1C − T1 ) 3.2. Analysis of the central part of the junction section As noted above, the expansion given by Eq. (26) diverges in the immediate neighborhood of the junction point C (i.e. for σ1 < ξ < σ2 , where (σ1 , σ2 ) → (0, 0) when δ → 0). In this region solution behavior can be examined by introducing new variables u, v, and τ as well as the reformulated expansion parameter μ: ξ = μu,

η = μ2 v,

t = μ2 τ,

μ3 = γδ.

(27)

γ(μu, μ2 v) . γ(0, 0)

(28)

Using these variables Eq. (20) becomes dv γ˜ (μu, μ2 v) = ; du u2 + v

γ˜ (μu, μ2 v) = 10

We now expand v(u) as follows: v(u) = v0 (u) + μv1 (u) + μ2 v2 (u) + μ3 v3 (u) + ... .

(29)

Insertion and ordering in powers of μ yields the following zeroth- to fourth-order equations: μ0 : μ1 : μ2 : μ3 :

μ4 :

1 dv0 = 2 , du u + v0 v1 dv1  + 2 = u˜ γξ , (u2 + v0 (u)) du u + v0 (u) v2 dv2 1  1 dv12  + 2 = v0 γ˜η + γ˜ξ2 u2 − , (u2 + v0 (u)) du u + v0 (u) 2 2 du v3 dv3 1   (3) + 2 = v1 γ˜η + γ˜ηξ uv0 + u3 γ˜ξ3 − (u2 + v0 (u)) du u + v0 (u) 6 d(v1 v2 ) , − du v4 dv4 1    + 2 = v2 γ˜η + γ˜ηξ uv1 + v02 γ˜ η2 + (u2 + v0 (u)) du u + v0 (u) 2 2 1 1 d(v1 v3 ) 1 dv2 (3) (4) + u2 v0 γ˜ ηξ2 + u4 γ˜ ξ4 − − . 2 24 du 2 du

(30) (31) (32)

(33)

(34)

The analysis is initiated by examining Eq. (30). By regarding u as a function of v0 , this becomes a Riccati equation: du = u 2 + v0 , dv0

(35)

which by the transformation u(v0 ) → −

1 df (v0 ) f (v0 ) dv0

(36)

assumes the form of Airy’s equation: d2 f (v0 ) + v0 f (v0 ) = 0. dv02 Following Watson [24], this has the solution   √ 2 3/2 2 3/2 f (v0 ) = v 0 C1 J−1/3 ( v0 ) + C2 J1/3 ( v0 ) . 3 3 11

(37)

(38)

Here Jp (z) is the order-p Bessel function of the first kind. Reverting to the dependent variable u(v0 ), we find for v0 ≥ 0: √ J−2/3 ( 23 v03/2 ) − cJ2/3 ( 32 v03/2 ) u(v0 ) = − v 0 , 3/2 3/2 cJ−1/3 ( 23 v0 ) + J1/3 ( 23 v0 )

(39)

where c = C1 /C2 . For v0 ≤ 0, an alternative representation of the solution to the Airy equation is useful, viz.  √ I−2/3 ( 23 (−v0 )3/2 ) − cI2/3 23 (−v0 )3/2   , u(v0 ) = − −v0 (40) cI−1/3 23 (−v0 )3/2 − I1/3 23 (−v0 )3/2 where Ip (z) is the order-p modified Bessel function of the first kind. Following [20], we note that for c = 1, u → +∞ as v0 → Ω0 , where Ω0 = 2.338 (known as the Haag special constant [22]) is the smallest value of v0 yielding 3/2 3/2 J−1/3 ( 23 v0 ) + J1/3 ( 23 v0 ) = 0. This specific solution, denoted the Haag special function (or the dividing solution) and starred, i.e. v0 , will play an important part in what follows and is shown in Fig. 3. The limiting behavior of the general solution u as v0 → −∞ is established by making use of Eq. (40) for u. To accomplish this, two useful relationships, valid for z → −∞, are   2 a1 (μ) 1/2 I−μ (z) ∼ Iμ (z) + sin(μπ)(π/2z) exp(z) 1 + + ... , (41) π z exp(z) Iμ (z) ∼ (2πz)1/2

 a1 (μ) + ... , 1− z



(42)

where heavily truncated versions of Hankel’s asymptotic expansions of the modified Bessel functions of the first and second kind have been employed. We are dealing with the cases μ = ±1/3, ±2/3 and hence a1 (±1/3) = −5/72,

a1 (±2/3) = 7/72.

(43)

When evaluating the limiting behavior of u as v → −∞ it is necessary to distinguish between two cases: When c = 1,  7 √ (1 − c) 1 − 48 √ (−v0 ))3/2 + ...  ∼ v0, lim u(v) = − v 0 (44) 5 3/2 v→−∞ (c − 1) 1 + 48 (−v0 )) + ... 12

and for c = 1,

 1+ lim u(v) = − v 0  v→−∞ 1+ √

7 (−v0 )3/2 48 5 (−v0 )3/2 48

√ + ... ∼ − v0. + ...

(45)

Thus, for c = 1 and v0 → −∞, all solutions u approach the upper branch of the parabola v0 + u2 = 0, whereas the only solution to approach its lower branch when v0 → −∞ is the one for c = 1, viz. the dividing solution v0 . 3.3. Behavior of the dividing solution The problem of determining the behavior of v0 (u) as u → +∞ is dealt with by reverting to Eq. (30), which after transformation by u → x−1 becomes dv  1 − 0 = ∼ 1 − v0 x2 . (46) dx 1 + v0 x2 Solving this equation is accomplished by the introduction of a power series v0 (x) =

∞

ak x k ,

(47)

k=0

which ultimately yields a0 = Ω0 , a1 = −1, a2 = 0, a3 = Ω0 /3, a4 = −1/4, a5 = 0, a6 = Ω0 /18. Here it should be noted that the limiting value Ω0 (= a0 ) is known from the analysis of the Riccati equation solutions in the previous subsection. Reverting to the original independent variable u it is found that Ω0 1 1 v0+ = Ω0 − + 3 − 4 + ...; u → +∞. (48) u 3u 4u The asymptotic behavior of the dividing solution v0 (u) as u → −∞ is found using an iterative procedure: (n−1)

dv0 du

=

1 u2 +

(n) v0

=⇒

(n) v0

=

1 − u2 (

(n−1)

dv0 du

(n−1) dv0

)

;

n = 1, 2, 3..., (49)

du

where v0 (u) = −u2 . Proceeding by using a truncated version of the formula for the summation of a geometric series at each step, the following result is ultimately obtained: v0− (u) ∼ −u2 −

1 1 5 − 4− ...; 2u 8u 32u7 13

u → −∞.

(50)

Once these limiting characteristics of v0 as u → ±∞ have been derived, it is possible to determine the analogous properties of the higher-order corrections governed by Eqs. (31)-(34). (For notational convenience, the dividing solution will henceforth be denoted v0 , i.e. stars will be dropped.) As u → +∞ the first-order equation to be resolved iteratively is represented by 

u˜ γξ dv1 (u)+ v1 (u)+ + . = 2 2 du u + v0 (u) (u2 + v0 (u))

(51)

Using the previously derived expression as well as the summation formula for a geometric series, this equation can be reformulated as     dv1+ v1+ 2Ω0 1 Ω0 2 3Ω20 1 Ω20 Ω0  + 4 1 − 2 + 3 + 4 + ... = γ˜ξ − + 4 + 5 − 6 + ... .(52) du u u u u u u3 u u 3u Here the dominant lowest-order (in terms of 1/u) balance can be expressed as (0)+

dv1 du



γ˜ξ = , u

(53)

whereafter an integration yields (0)+

v1



(u) = γ˜ξ ln u + Ω1 .

(54)

The integration constant Ω1 characterizes the specific property of v1+ distinguishing this function in the set of all solutions to Eq. (31). It is now assumed that (1)+ v1 (u)

(1)+



= γ˜ξ ln u + Ω1 + f1 (u)

dv =⇒ 1 du



γξ df1 = + . u du

(55)

Following the same path as above ultimately yields 



γ ξ Ω0 γ ξ Ω0 df1 = − 3 =⇒ f1 (u) = + f2 (u), du u 2u2

(56)

and thus 

(1)+ v1 (u)

γ ξ Ω0 = γ˜ξ ln u + Ω1 + + f2 (u). 2u2 

14

(57)

The change of variables given by Eq. (27), together with our choice to restrict the present analysis to the first-order expansion in terms of δ, makes the term 

γξ Ω 0 2u2

+ f2 (u) redundant, and thus (1)+

v1



(u) = γ˜ξ ln u + Ω1 + ... ,

(58)

whereafter the process can be repeated to resolve Eq. (52) to any desired level of approximation. The constant Ω1is obtained by numerical evaluation  (1)+ of the quantity lim v1 (u) − γ˜ξ ln u , yielding Ω1 = 3.950. u→+∞

Subjecting Eqs. (32)-(34) to an analogous treatment, and using the same argument as above for restricting the analysis to the power μ3 (equivalent to δ), yields the following second- to fourth-order results: 

v2+ (u)

γξ 2

u + ..., 2 (3) γξ 3 2 + v3 (u) = u + ..., 12 (4) γξ 4 3 + v4 (u) = u + ..., 72 =

(59) (60) (61)

which ultimately yields 

1 γξ ln δγ − T2 (T1 , δ) = T2C − δ 2/3 γ 2/3 Γ − 31 δγ˜

 −δγ

−M (T1C −T1 )



+ γ˜ξ ln

T1C −T1 M

+ Ω1 +

γ ˜



ξ2

(T1C −T1 ) 2M

+

γ ˜

(3)

ξ3

(T1C −T1 )2 12M 2

+

γ ˜

(4)

ξ4

(T1C −T1 )3



72M 3

+ ... . (62)

In the limit of u → −∞ a similar iterative procedure can be applied when dealing with Eqs (31)-(34). Making use of Eq. (50), the higher-order corrections are in this limit found to be 

γ˜ξ = − + ..., 2   2˜ γ ˜ξ η −γ − u + ..., v2 (u) = 4  (3) 6˜ γηξ − γ˜ξ 2 − v3 (u) = u + ..., 12 v1− (u)

15

(63) (64) (65)

(3)

v4− (u)

=



(4)

12(˜ γηξ2 − γ˜η2 ) − γ˜ξ4 48

u3 + ... .

(66)

Once these results for u → ±∞ have been established, it is possible to determine the functions vn (u); n = 1, ..., 4 in their entirety. This was accomplished using a five-step Kutta-Merson numerical scheme and employing the asymptotically valid results as u → −∞ for initialization purposes. These results, to be used when constructing the approximated limit cycle, were validated on the basis of the results above for vn+ (u); n = 1, ..., 4 as u → +∞ and are shown in Fig. 3, which also includes the Haag special function v0 (u). 3.4. Analysis of the end of the junction section The expansion given by Eq. (62) is not valid over the entire segment b1 ≤ ξ ≤ b2 . We thus revert to Eq. (20), which is reformulated using the expansion parameter μ: dη γ˜ (ξ, η) = μ3 2 . dξ ξ +η

(67)

In view of the results of the previous subsection, the following expansion of η can be assumed: 1 η = μ2 Γ0 (ξ) + μ3 ln Γ1 (ξ) + μ3 Γ2 (ξ) + ... . μ

(68)

Insertion and ordering in powers of μ yields the following second- and higherorder problems: dΓ0 = 0 =⇒ Γ0 = Ω0 , dξ dΓ1 1  : = 0 =⇒ Γ1 = γ˜ξ , μ3 ln μ dξ  ξ γ˜ (x, 0) γ˜ (ξ, 0) dΓ2 3 = =⇒ Γ = dx + Ω1 . μ : 2 dξ ξ2 x2  μ2 :

(69) (70) (71)

In these solutions the constants of integration have been chosen so as to ensure matching with Eq. (62). The integral in Eq. (71) is obviously divergent as → 0+ , but can be dealt with using an approach due to Hadamard [25].

16

It is thus reformulated as



ξ

P (ξ) =  

ξ

1 x2



(m) 1

γ˜ξm (0, 0) m x γ˜ (x, 0) − m! m=0







ξ

dx + 

γ˜ (x, 0) dx = x2

(m) 1

γ˜ξm (0, 0) dx. m!x2−m m=0

An integration of the first r.h.s term yields (m+1) (m+1) (m)  ξ ∞ ∞ ∞

γ˜ξm+1 (0, 0) m γ˜ξm+1 (0, 0) m γ˜ξm (0, 0) m−2 x ξ −

. dx = m! m(m + 1)! m(m + 1)!  m=2 m=1 m=1

(72)

(73)

The second of the r.h.s. terms does not cause any problems since it approaches zero as → 0+ . Integrating the remaining term of P (ξ) yields     (m)  ξ 1 γ˜ξm (0, 0) 1 1   − γ˜ξ (0, 0) ln + − + γ˜ξ (0, 0) ln ξ , (74) dx = 2−m

ξ  m=0 m!x where the - dependent term evidently diverges as → 0+ and thus can be neglected, this since we only make use of the finite terms of the asymptotic expansion of P (ξ), viz. 

(3)

(4)

γ˜ξ3 2 γ˜ξ4 3 γ˜ξ2 1  Pf in (ξ) = − + γ˜ξ ln ξ + ξ+ ξ + ξ ... . ξ 2 12 72

(75)

The desired expansion to first order in δ thus assumes the following form: T2 (T1 , δ) = T2C − δ 2/3 γ 2/3 Ω0 − 1  1 − γ˜ξ δγ ln − δγ(Ω1 + Pf in (ξ)) + ... . 3 δγ

(76)

This result is evaluated at T1 = b2 and used for initializing the subsequent part of the limit cycle characterized by fast phase-point motion. 4. Rapid motion towards the drop region A new expansion is required for a proper description of this rapid-motion segment of the limit cycle, which, as noted above, is taken to commence at T1 = b2 and is governed by the equation 1 dT2 = Hν (T1 , T2 ); δ dT1 17

ν = 1, 2.

(77)

That ν assumes both of its values reflects that the rapid motion takes place in both phase-plane sectors T22 > T12 and T22 < T12 . To proceed, the solution is expanded in powers of δ 1/3 and ln(1/δ): (0)

(2/3)

T2 (T1 , δ) = T2 (T1 ) + δ 2/3 T2

1 (2/3) (1) (T1 ) + δ ln T2 (T1 ) + δT2 (T1 ) + ..., (78) δ

where experience has proved that including a term of order δ 1/3 is redundant. Insertion and ordering yields the following sequence of problems: dT2 (0) = 0, dT1 dT2 (0) 2/3 δ : = 0, dT1 1 dT2 (2/3) δ ln : = 0, δ dT1 dT2 (1) (0) δ1 : = Hν (T1 , T2 (T1 )); dT1 δ0 :

(79) (80) (81) ν = 1, 2.

(82)

This system of differential equations is integrated and taken to be constrained by the evaluation of T2 (b2 , δ) based on the results of the previous section, which yields δ 0 : T2 (0) (T1 ) = T2C (= T2D ), δ 2/3 : T2 (2/3) (T1 ) = −γ 2/3 Ω0 , 1 1  δ ln : T2 = − γ˜ γ, δ 3 ξ  T1 1  1 (1) 1 δ : T2 (T1 ) = − γ˜ξ γ ln − Ω1 γ − Pf in (b2 )γ + Hν (x, T2C )dx; 3 γ b2

(83) (84) (85) ν = 1, 2. (86)

The expansion of T2 (T1 , δ) along this rapid-motion segment of the limit cycle is thus given by 

T1

T2 (T1 , δ) = T2 (b2 ) + δ b2

H1 (z, T2C )dz + ... = T2D + Ξ(T1 , δ),

where T2 (b2 ) is evaluated using Eq. (76). 18

(87)

 T1

Hν (x, T2C )dx; ν = 1, 2 to a very good approximation  T1 could be simplified to b2 H1 (x, T2C )dx is demonstrated in Appendix B. This rapid-motion segment is hence to first order given by  T1 −z 3 + k1 z 2 + k2 z + k3 dz, (88) T2 (T1 , δ) = T2 (b2 ) + δ z 3 − k1 z 2 − k4 z + k5 b2 That the integral

b2

where k1 = T2C , k2 = κ/Ra + k12 , k3 = 1/[(1 − θ)Ra] − (κ + 1)k1 /Ra − k13 , k4 = (κ + 1)/Ra + k12 , k5 = θ/[(1 − θ)Ra] + κk1 /Ra + k13 . The expression for T2 (T1 , δ), evaluated at the end of the rapid-motion segment, will be essential for initiating the analysis of solution behavior in the region adjacent to the drop point (T1D , T2D ), the location of which is derived in Appendix A. From the definition of H2 (T1 , T2 ) given by Eq. (3) it is recognized that for T1 → T1D , the first-order correction becomes unbounded. Hence the validity of Eq. (88) is only taken to extend to T1 = b3 . This point, located within the region of attraction of F (T22 < T12 ) = 0 and to be determined when representing the approximated limit cycle, is where the analysis of solution behavior in the drop region is to be initiated. 5. The drop region adjacent to D In this region, special measures must be employed. The analysis  thus proceeds from an ”auxiliary” drop point D† , viz. b3 , T2D + Ξ(b3 , δ) in the vicinity of D. A local-variable system (p, q) is introduced by the following transformations: T1 = T1D + p(q, δ),

T2 = T2D + Ξ(b3 , δ) + δβD ln q,

t = δθ,

(89)

where βD = −G(T1D , T2D )/FT 1 (T1D , T2D ). Substitution into δ

dT1 = H2−1 (T1 , T2 ) (= M2 (T1 , T2 )) , dT2

(90)

yields the following equation: dp βD = M2 (T1D + p, T2D + Ξ(b3 , δ) + δβD ln q), dq q 19

(91)

constrained by p(1, δ) = b3 − T1D << 1. We now seek the solution to this problem in the form of an expansion in powers of δ 1/3 and ln(1/δ) : 1 p = p0,0 + δ 2/3 p2,0 + δ ln p3,1 + δp3,0 + ... , δ

(92)

1 Ξ = δ 2/3 Ξ2,0 + δ ln Ξ3,1 + δΞ3,0 + ... , δ

(93)

together with

where Ξ2,0 = −γ 2/3 Ω0 , 1  γ, Ξ3,1 = − γ˜ 3 ξ  b3 G(z, T2C ) 1  1 dz. Ξ3,0 (q3 ) = − γ˜ξ γ ln − γΩ1 − Pf in (b2 )γ + C 3 γ b2 F (z, T2 )

(94) (95) (96)

Here T1 = b3 is, as mentioned in the previous section, the limit where the fast motion is taken to end. As for T2 (T1 , δ) over the rapid-motion segment, it was a posteriori found superfluous to include a term of order δ 1/3 in this expansion of p. Insertion and ordering yields the following sequence of problems: βD dp0,0 = M2 (T1D + p0,0 (q), T2D ), dq q dp2,0 βD   D − M2T1 T1 + p0,0 (q), T2D p2,0 (q) = δ 2/3 : dq q βD Ξ2,0  M2T2 (T1D + p0,0 (q), T2D ), = q dp3,1 βD  1 : − M (T D + p0,0 (q), T2D )p3,1 (q) = δ ln δ dq q 2T1 1 βD Ξ3,1  M2T2 (T1D + p0,0 (q), T2D ), = q dp3,0 βD  − M2T1 (T1D + p0,0 (q), T2D )p3,0 (q) = δ1 : dq q βD (Ξ3,0 (b3 ) + βD ln q)  = M2T2 (T1D + p0,0 (q), T2D ), q δ0

:

20

(97)

(98)

(99)

(100)

subject to the conditions p0,0 (1) = b3 − T1D ,

p2,0 (1) = p3,1 (1) = p3,0 (1) = 0.

(101)

Making use of the formal definition of M2 (T1 , T2 ), an integration of the zerothorder problem yields the desired solution    p0,0 G(T1D + x, T2D ) 1 q(p0,0 ) = exp dx , (102) βD b3 −T1D F (T1D + x, T2D ) which in turn can be inverted to provide p0,0 (b). Once this is known, the remaining functions can be successively derived:  D  q  T1 + p0,0 (x), T2D Ξ2,0 M2T 2 p2,0 (b) = βD Π(q) dx, (103) xΠ(x) 1  p3,1 (b) = βD Π(q)  p3,0 (b) = βD Π(q)

q 1

q 1

 D  D T Ξ3,1 M2T + p (x), T 0,0 1 2 2 dx, xΠ(x)

(104)

 D  D (Ξ3,0 (b3 ) + βD ln x)M2T + p (x), T T 0,0 1 2 2 dx, (105) xΠ(x)

where  Π(q) = exp 1

q

 βD M2T (T1D + p0,0 (x), T2D ) 1 dx. x

(106)

The higher-order approximation of this trajectory segment is hence given by the following expansions: T1 (q, δ) = T1D + p0,0 (q) + δ 2/3 p2,0 (q) + δ ln 1δ p3,1 (q) + δp3,0 (q) + ... ,(107) T2 (q, δ) = T2D + δ 2/3 Ξ2,0 + 13 δ ln 1δ Ξ3,1 + δ (Ξ3,0 (b3 ) + βD ln q) + ... .(108) The evaluation of these expressions is pursued for increasing values of q until a smooth merger with the slow-motion segment near F (T22 < T12 ) = 0 has been accomplished.

21

6. Slow motion in the sector T22 < T12 In contrast to how the analysis of the slow-motion segment in the proximity of F (T22 > T12 ) = 0 was conducted in Section 2, T1 will here be regarded as a function of T2 and use will thus be made of Eq. (1) divided by Eq. (2): δ

dT1 = H2−1 (T1 , T2 ). dT2

To resolve the problem, the following power series is introduced: (0)

(1)

T1 (T2 , δ) = T1 (T2 ) + δT1 (T2 ) + ... ,

(109)

whereby Eq. (109) can be reformulated as (0)

(1)

dT AT + δBT + ... dT . δ( 1 + δ 1 + ...) = dT2 dT2 AN + δBN + ...

(110)

The analysis pertains to the sector T22 < T12 , and thus (0) 2

(0)

(0) 2

(0)

(0)

(0)

θ AT = −Ra(T22 − T1 )(T2 − T1 ) + ( 1−θ − T1 ) + κ(T2 − T1 ), (111)   (1) (0) 2 (0) (112) BT = T1 Ra(T22 − 3T1 + 2T2 T1 ) − (κ + 1) , (0)

1 )(T2 − T1 ) + ( 1−θ − T2 ) − κ(T2 − T1 ),   (1) (0) 2 (0) BN = −T1 Ra(T22 − 3T1 + 2T2 T1 ) − κ .

AN = Ra(T22 − T1

(113) (114)

Reordering the problem in powers of δ yields 0

AT = 0,

δ :

(0)

dT AN 1 = BT . dT2

1

δ :

(115)

Differentiating the lowest-order result with respect to T2 yields (0) 2

(0)

(0)

Ra(T1 + 2T1 T2 − 3T22 ) + κ dT1 , = (0) 2 (0) dT2 Ra(3T1 − T22 − 2T1 T2 ) + κ + 1

(116)

whereby the first-order result becomes (0)

(1) T1 (T2 )

=

dT1 dT2 (0)

AN (0) 2

Ra(2T1 T2 + T22 − 3T1

) − (κ + 1)

.

(117)

As to be discussed further in Section 10, the expansion based on this first-order correction term converges extremely rapidly. 22

7. Solution behavior subsequent to reaching the cusp In its slow motion near F (T22 > T12 ) = 0, the zeroth-order solution loses stability as it ultimately reaches the cusp A where the two branches of F (T1 , T2 ) = 0 merge on the antisymmetric diagonal, cf. Appendix A. (This instability is the subject of a forthcoming discussion in Section 10.) Similar behavior is found when the first-order correction term, cf. Eq. (117), is taken into account. In this case, however, instability sets in at the modified ∗ cusp point A∗ . This is found by determining T2A from the requirement that A∗ also be located on the antisymmetric diagonal, viz. ∗

∗

T2A = −T1 (T2A , δ).

(118)

Reverting to the first-order relationship in Eq. (117) above, it is recognized (0) that AN ∼ dT2 /dt and dT1 /dT2 are positive and negative definite, respectively, cf. the phase-plane representation of the limit cycle. However, the denominator of the ratio constituting the correction term is negative defi(1) nite, and thus T1 (T2 ) > 0. This implies that A∗ is shifted towards the origin along the antisymmetric diagonal compared to the location of the lowest-order cusp A. The initial rapid motion from A∗ in the sector T22 > T12 is dealt with by reformulating the governing equations as 1 dT2 = H1 (T1 , T2 ), δ dT1

(119)

and assuming that (0)

(1)

T2 (T1 , δ) = T2 + δT2 (T1 ) + ... .

(120)

Insertion and ordering in powers of δ yields the results already reported in Section 2, cf. Eq. (7). However, since KN = 0 during this rapid phase, dT

(0)

the lowest-order result assumes the form dT21 = 0, with concomitant effects on the first-order equation. Both relationships can trivially be integrated to yield (0)

δ 0 : T2 δ1

∗

= T2A ,  (1) : T2 (T1 ) =

(121) T1 ∗

T1A

3



2





−z + k1 z + k2 z + k3    dz, z 3 − k1 z 2 − k4 z + k5 23

(122)

where 

∗

k1 = T2A , 



2

k2 = κ/Ra + k1 , 2

k4 = (κ + 1)/Ra + k1 ,







3

k3 = 1/[(1 − θ)Ra] − (κ + 1)k1 /Ra − k1 , 

3

k5 = θ/[(1 − θ)Ra] + κk1 /Ra + k1 .

Since the drop point (T1B , T2B ) is known, cf. Appendix A, it is recognized that as T1 approaches T1B , this correction term of order delta becomes unbounded. Hence Eq. (120) is only taken to be valid up to T1 = b4 , located within the domain of attraction of F (T22 > T12 ) = 0. Thus, as recognized in Section 5, when the solution approaches the drop region a situation arises which requires to be dealt with in a less straightforward manner. 8. The drop region adjacent to B The analysis in this region is initiated by focusing on the trajectory segment derived in the preceding section, evaluated at an ”auxiliary” drop point   † (1) † B B located a small finite distance from B with coordinates b4 , T2 + δT2 (T1B ) . We use basically the same procedure as in Section 5, but here applied to the equation δ

dT1 = H1 (T1 , T2 ) (= M1 (T1 , T2 )) . dT2

Local ”drop-region” variables (z, w) are introduced, viz.   (1) B B B† T1 = T1 − z, T2 = T2 + δ T2 (T1 ) + βB ln w ,

(123)

(124)

where βB = −G(T1B , T2B )/FT 1 (T1B , T2B ). Substitution of these variables into the pertinent form of Eq. (3) yields   dz βB † (1) = − M1 T1B − z, T2B + δ(T2 (T1B ) + βB ln w) , dw w

(125)

†

subject to the constraint z(1, δ) = T1B − T1B << 1. We seek a solution in the form of a power series z(w, δ) = z (0) (w) + δz (1) (w) + ... . Insertion and ordering in powers of δ yields the following zeroth- and first-order problems:

δ0

:

βB dz (0) = − M1 (T1B − z (0) (w), T2B ), dw w 24

(126)

δ

1

dz (1) βB   B − M1T1 ( T1 − z (0) (w), T2B z (1) (w) = : dw  w  † (1) βB T2 (T1B ) + βB ln w  B  = − M1T T1 − z (0) (w), T2B , 2 w

(127)

†

where z (0) (1, δ) = T1B − T1B and z (1) (1, δ) = 0. The zeroth-order problem has the solution    z(0) B B G(T1 − x, T2 ) 1 w(z (0) ) = exp − dx , (128) † βB T1B −T1B F (T1B − x, T2B ) inversion of which yields z (0) . Use of the integrating factor  Π(w) = exp 1

w

 B  (0) B T βB M1T − z (x), T 1 2 1 dx x

(129)

leads to the following first-order solution:    w T (1) (T B † ) + βB ln w M  T B − z (0) (x), T B 2 1 2 1 1T2 dx. z (1) (w) = βB Π(w) xΠ(x) 1 (130) Consequently the results valid to first order for this drop-region segment of the approximated limit cycle are T1 (w, δ) = T1B − z (0) (w) − δz (1) (w) + ... ,   † (1) T2 (w, δ) = T2B + δ T2 (T1B ) + βB ln w + ... .

(131) (132)

As in Section 5, these expressions are evaluated for an increasing w until a smooth merger with the slow-motion segment has been achieved. Since the subsequent slow motion in the sector T22 > T12 towards the junction point C already has been dealt with in Section 2, this set of results brings the present analysis of the limit cycle associated with Eqs. (1) and (2) to en end.

25

9. Results Evaluations based on the procedures described above have been carried out to determine the asymptotically valid corrections to the discontinuous zeroth-order solution of Eqs. (1) and (2). To facilitate comparisons with previous results [17], the choice of parameters was Ra = 10, θ = −13/14, and κ = 0.05. Two values of the expansion parameter were employed, δ = 0.1 and 0.05. From a procedural standpoint it should be noted that, although the general thrust of the present investigation is analytical, the integrals occurring in certain of the formulae in Sections 3-8 have been evaluated numerically. Furthermore, the ”auxiliary” drop points D† and B † (cf. Sections 5 and to initiate the drop regions are chosen such that   8) introduced  D − D†  << 1 and B − B †  << 1. Since the analysis of the limit cycle has proceeded segment by segment making use of the results in Appendix C, it should also be emphasized that the associated matching in general constitutes an integral part of the analytical procedure. In the sector T22 > T12 the merger between the initial part of the junction section and the approximated slow-motion segment based on a correction of order δ was not entirely satisfactory. There is a possibilty of ameliorating this shortcoming by extending the slow-phase analysis to order δ 2 on the basis of the second-order relationship from Eq. (7) and the explicit versions of LT and MN given in Eq. (6). These formulae as well as previously derived lower-order results yield δ2 :

(2) T2 (T1 )

(1) 2

=−

T2

(0)

(1)

Ra(3T2 − T1 ) + T2 − LN ( (0)

(0) 2

Ra(2T1 T2 + T2

(1)

dT2 dT1

+ 1)

− 3T12 ) + κ + 1

, (133)

where (1)

dT2 1+θ (0) (0) (0) 2 (3T1 − T2 ) − 3T12 − 6T1 T2 + T2 + = [Ra{2 dT1 1−θ (0) 1+θ dT (0) (0) 2 (0) (T1 + T2 ))} − + 2 (5T12 + 2T1 T2 + T2 − 2 dT1 1−θ (0)   −2  dT (0) 2 (0) +κ+1 . −(κ + 1)(1 + 2 )] Ra −3T12 + T2 + 2T1 T2 dT1

(134)

As a contrast, it deserves note that the behavior of the initial part of the junction section is almost impervious to an extension of the order of analysis. 26

Phase-plane representations of the asymptotically valid approximations of the relaxation oscillation are shown as solid lines in Figs. 4 and 5 for δ = 0.05 and 0.1, respectively. The graphs also include the corresponding numerical (dash-dotted) and discontinuous lowest-order (dashed) results previously reported in [17]. An overall judgement of the asymptotically valid corrections derived in the present study is that they serve a useful purpose. As expected, and confirmed by a comparison between Figs. 4 and 5, their fidelity to the numerical results decreases as the expansion parameter δ grows larger, but nevertheless these solutions must be regarded as considerable improvements upon the previously derived zeroth-order results in [17] and [21]. 10. Physical mechanisms underlying the relaxation oscillation An alternative overview of the behavior of the system is provided by Fig. 6 pertaining to the ”standard” case of Ra = 10, θ = −13/14, κ = 0.05, and δ = 0.1. This shows the numerically calculated time evolution of T1 , T2 , and T2 + δT1 , the latter quantity being the total heat content of the system. Fig. 7 analogously shows the behavior of the heat fluxes acting on and within the system. These are the advective flux through the connecting tubes, the conductive fluxes through the outer walls of the vessels (forced by the prescribed external temperatures T1ex and T2ex ), and the conductive flux through the internal common internal wall of the system. From these diagrams, as well as from the phase-plane representations of the limit cycle, it is recognized that on the slow-motion segment in the sector T22 < T12 , T1 and T2 are located comparatively close to the antisymmetric diagonal. During this slow-motion phase, the dynamics of the system are thus dominated by the conductive heat fluxes through the outer walls of the vessels since only a small advective heat exchange takes place. To estimate (1) the magnitude of the correction term T1 (T2 ) to be applied to this segment of the limit cycle, Eq. (117) is reformulated as (1)

T1 (T2 ) =

  (0) (0) (0) (0) ((1 + θ)/(1 − θ) − (T2 + T1 ))(Ra (T2 − T1 )2 + 2(T2 − T1 )(T2 + T1 ) − κ) .  2 (0) 2 (0) (0) Ra(T2 − T1 ) − 2(T2 − T1 )(T2 + T1 ) + κ + 1 (135) 27

(0)

(0)

Since T2 − T1 > 0 and T2 + T1 < 0 during this slow phase near F (T22 < T12 ) = 0 , a ”reduced” version of the equation above can be used to establish an upper bound on the magnitude of the correction term: |

(1) T1 (T2 )

(0)

|<

(1 + θ)/(1 − θ) − (T2 + T1 ) (0)

(Ra(T2 − T1 )2

.

(136)

When evaluating this expression for θ = −13/14 and Ra ranging between 10 (0) and 250, a conservative estimate is that T1A < T1 < T1D and T2D < T2 < T1A . It is found that the correction term in principle can be disregarded when representing the limit cycle, this since its maximal absolute values for small and large Ra are on the order of 10−2 and 10−3 , respectively. The prescribed thermal forcing of the system is such that | T1ex /T2ex |< 1, and thus the slow motion in the sector T22 < T12 is inexorably directed towards the antisymmetric diagonal, which ultimately is reached at the cusp   (T1A , −T1A ). Here no advective exchange between the vessels takes place and     although F (T1A , −T1A ) = 0, this state is not stationary since G(T1A , −T1A ) = (1 + θ)/(1 − θ). This small conductive heat flux into the larger vessel can be interpreted as resulting in a temperature perturbation ΔT2 > 0, which, however, is not destabilizing as shown by a linear stability analysis at the cusp   (T1A , −T1A ) based on Eq. (2). A similar analysis of Eq. (1) for (T1 +ΔT1 , T2 ), where ΔT1 is taken to be positive and assumed to be proportional to eλt , yields 

4Ra(T1A )2 − (κ + 1) λ= . δ

(137)

Since λ > 0 for the parameter values under consideration in this study, the solution at the cusp is unstable for small positive perturbations of T1 , whereas it is unconditionally stable for ΔT1 < 0. The momentary equilibrium on the antisymmetric diagonal precludes any advective fluxes, but any small perturbation of T1 will alter the conductive heat fluxes affecting the smaller vessel and also induce advective heat fluxes tending to increase T1 and reduce T2 . A negative perturbation of T1 will thus be counteracted, whereas a positive perturbation will reinforced, resulting in a loss of stability. When instability sets in at the cusp, T1 increases due to the advective heat flux from the larger vessel to the smaller, which dominates over the conductive heat loss from the latter vessel. Initially T2 also shows a small   increase since G(T1A , −T1A ) > 0, where the externally forced heat flux to the 28

larger vessel briefly dominates over the advective flux of colder water to this vessel. This initial stage comes to an end when the advective and conductive heat fluxes affecting the larger vessel balance one another, whereafter the former flux prevails. The smaller vessel is now flushed by water from the larger; this rapid-motion phase is initially dominated by the advective heat exchange, which, although | T22 − T12 | only is moderately large near T1 = −T2 , nevertheless is important since T2 − T1 assumes a significant magnitude here. This increases the smaller-vessel temperature T1 despite a growing conductive heat loss from this vessel through its outer wall. As the slowmotion segment in the sector T22 > T12 is approached, the imbalance between these two dominating heat fluxes affecting T1 gradually decreases and the rapid-motion of the phase-point slows down as dT1 /dt diminishes towards its turning-point value of zero on F (T22 > T12 ) = 0. At this stage the conductive heat loss from the smaller vessel is of a significant magnitude, since here T1 −T1ex assumes its maximum value during the entire limit cycle. At the same time, the amount of heat gained by the larger vessel from thermal conduction through its outer wall is comparatively modest, since here T2ex − T2 assumes values not far removed from its minimum value at the cusp. In conjunction with the advective flux of colder fluid from the smaller vessel this has the consequence that T2 decreases, hereby giving rise to a smooth merger with the slow-motion segment. As recognized from Fig. 7, the physical mechanisms underlying the behavior of the solution during its subsequent slow motion in the sector T22 > T12 are of a somewhat more subtle nature than was the case during the analogous phase in the other sector, which, as seen above, was dominated by the externally forced conductive heat fluxes. As seen in Figs. 4 and 5, when T22 > T12 the slow-motion segment is located at a considerable distance from the symmetric diagonal T1 = T2 . This indicates that here the advective heat fluxes a priori can be expected to play a significant role. On this segment, T1 decreases since the conductive heat loss from the smaller vessel through its outer wall (T1 >> T1ex ) exceeds the advective flux of warmer water from the larger vessel. Since δ << 1, immediately after the ”jump” T2ex − T2 is comparatively small, and thus the externally forced heat flux to the larger vessel is overshadowed by the advective flux of colder water from the smaller vessel, which makes also T2 decrease. As both T1 and T2 continue to decrease, the advective heat exchange between the vessels becomes smaller at the same time as the conductive fluxes affecting the larger and smaller vessel vessel wax and ∗ ∗ wane, respectively. At (T1C , T2C ) in the neighborhood of the lowest-order 29

junction point C, the system ultimately reaches a state dT2 /dt = 0 when the conductive heat fluxes affecting the larger vessel momentarily balance the effects of the advection of colder water from the smaller vessel. (This situation corresponds to the minimum of T2 in the phase-plane, since over the greater part of the limit cycle a decreasing T1 can be interpreted as proxy for t.) ∗ ∗ At (T1C , T2C ) Eqs. (1) and (2) thus simplify to δ

dT1 1+θ = − T1 − T 2 , dt 1−θ ∗

(138)

and by the substitution T1 = T1C ± ΔT1 , ΔT1 > 0 (it being assumed that ∗ ∗ ΔT1 ∝ eλt ), a stability analysis yields that λ = ∓ 1δ . At (T1C , T2C ) the solution thus becomes unstable for small negative perturbations of T1 . Here the mechanism underlying the loss of stability is somewhat less straightforward  than was the case at the cusp. It hinges on the fact that T1C < 0, and thus a positive perturbation of T1 yields a smaller absolute value of this temperature, which leads to an increased advective heat flux to the smaller vessel. Concurrent with this heat gain, however, an increased conductive heat loss through the outer wall of this vessel takes place, which is of such a magnitude  so as to restore T1 towards its equilibrium at T1C . A negative perturbation of T1 , on the other hand, decreases the advective flux of warmer water to the smaller vessel as well as its externally forced conductive heat loss. Since the effects of the altered advective flux dominate, T1 continues to decrease, i.e. instability sets in and the transition phase commences. During the ”jump” the same mechanism as that underlying the loss of stability continues to act, with the advective exchange gradually decreasing until it ceases at T1 = −T2 . Once this antisymmetric diagonal has been transgressed, the advective heat gain experienced by the smaller vessel begins to grow until it balances the conductive heat loss from this vessel, which concurrently has become smaller as T1 has decreased towards T1ex . For evident reasons this description of the transition phase has its focus on the behavior of T1 , but it is also noteworthy that T2ex − T2 remains comparatively large as this rapid phase comes to an end. At this stage the highly significant conductive heat flux to the larger vessel is of crucial importance for the behavior of the system, this since it has the consequence that the slow-motion phase in the sector T22 < T12 is initiated in a smooth fashion. Note finally that for conceptual clarity the only marginal (κ = 0.05 ) heat flux through the common wall separating the vessels has not been dealt with explicitly above. 30

11. Behavior of the system for Ra >> 1 The low value of the Rayleigh-type number previously dealt with, Ra = 10, corresponds to a weak thermal forcing of the system. Given that the physical properties of the working fluid can be regarded as immutably fixed, a stronger forcing, viz. a larger value of Ra, can most conveniently be accomplished by increasing the difference between the external-forcing temperatures and/or enlarging the diameter of the tubes connecting the two vessels. For Ra = 250, θ = −13/14, κ = 0.05, and δ = 0.1 the lowest-order as well as the asymptotic-analysis results show a remarkably enhanced adherence to the numerically calculated limit cycle, cf. Fig. 8. That the lowest-order and numerical results for strong forcing, viz. Ra >> 1, can be expected to coincide during the slow phases is recognized by dividing Eq. (1) with Ra. The left-hand member of the resulting equation is now constituted by the time-derivative of T1 preceded by the exceedingly small quantity δ/Ra and can thus be neglected. The remainder of the equation constrains T1 in terms of T2 and is evidently equivalent to the slow phases of motion as given by the lowest-order relationships KN = 0 and AT = 0. The overall problem remains well-posed since the evolution in time of the system is now governed by a first-order O.D.E. given by Eq. (2), where T1 has been expressed in terms of T2 based on the degenerate version of Eq. (1). As already noted in Section 6, for T22 < T12 the agreement between the numerical and the lowest-order solutions is very good even for small values of Ra. The improved global representation of the limit cycle for Ra >> 1 is thus primarily due to solution behavior in the sector T22 > T12 . A closer examination of a slightly reformulated version of the first-order correction term for the slow motion here reveals why this is the case: (1)

T2 (T1 ) =

(0)

(1 + θ)/(1 − θ) − (T2 + T1 )  . (0) (0) (0) κ + 1 − Ra (T2 − T1 )2 − 2(T2 − T1 )(T2 + T1 )

(139)

Large values of Ra diminish the overall magnitude of this term, cf. Fig. (1) 9 showing T2 (T1 ) versus T1 − T1C for different values of Ra. In this figure another feature of the first-order correction is also visible, namely that (1) T2 (T1 ) grows beyond bounds as T2 → T2C and T1 → T1C , cf. Section 2 . This holds true whatever the magnitude of Ra; the underlying mechanism is

31

conveniently illustrated by noting that for Ra >> 1, (1)

| δT2 (T1 ) |∝

1 1 δ   . (0) (0) Ra T − T T + 3T 1 1 2 2

(140)

T1C as a function of T2C is given by Eq. (A.3) in Appendix A, and for Ra >> 1 TC this relationship reduces to T1C ∼ − 32 , hereby explaining the singularity of (1) T2 (T1 ) as T1 approaches T1C . It is also recognized from Fig. 9 that for small Ra, this singular behavior starts to manifest itself at a considerable distance (1) from T1C , whereas for Ra >> 1, T2 (T1 ), remains minuscule until it ”spikes” in the immediate neighborhood of T1C . This feature has the beneficial consequence that for Ra >> 1, the somewhat subtle junction-section analysis (cf. Section 3) only needs to be implemented for a very minor segment of the limit cycle. In view of these considerations it is recognized that for large values of Ra, the asymptotically valid representation of the limit cycle should also be useful for larger values of δ than these hitherto dealt with, i.e. δ = 0.1 and δ = 0.05. This expectation is borne out by Fig. 10 showing the results for Ra = 250, θ = −13/14, κ = 0.05, and δ = 0.2. An interesting feature of Figs. 8 and 10 is that for these large values of Ra, the slow motion in the sector T22 > T12 initially appears to take place along a trajectory parallel to the symmetric diagonal T1 = T2 . The background to this behavior can be understood by examining the condition Kn = 0, the lowest-order representation of this segment of the limit cycle, since (as noted above) the first- and second- order corrections are of negligible magnitude for Ra >> 1. This equation can be reformulated as (0) 2

(T2

(0)

− T12 )(T2 − T1 ) + (

θ (0) − T1 ) − κ(T2 − T1 ) = 0, 1−θ (0)

where the parameter = 1/Ra is assumed to be small. T2 as (0)

T2

(0,0)

= T2

(141)

is thus expanded

+ 1/2 P (T1 ) + Q(T1 )... .

(142)

Insertion of this expression and ordering in powers of ultimately yields the following results: (0,0)2

0 : (T2

(0,0)

− T12 )(T2

− T1 ) = 0 32

=⇒

(0,0)

T2

= T1 ,

(143)



1

3/2

1/2  1 θ/(1 − θ) : P (T1 ) = √ 1 − , T1 2 (2κ + 1)T1 − θ/(1 − θ) : Q(T1 ) = − . 8T12

(144) (145)

The two correction terms suffice to include the lowest-order effects of the (0) parameters θ and κ in the series representation of T2 . Making use of this expansion it is found that for large values of T1 , (0)

T2 − T1 ∼

1/2 1 √ . = 21/2 2Ra

(146)

We furthermore note that (0)

dT2 θ =1+ √ 2 dT1 2 2T1 (1 − θ)



dT

θ 1+ 2T1 (1 − θ)

−1/2



1/2

+

2κ + 1 θ/(θ − 1) − 8T12 4T13

(0)

Hence, given Ra >> 1, dT21 → 1 for large values of T1 . This result is consonant with the slow-motion in the phase-plane sector T22 > T12 initially being parallel to the symmetric diagonal. Since the quantities P (T1 ) and Q(T1 ) furthermore increase in magnitude as T1 decreases, the lowest-order (0) 1 slow-motion trajectory begins to deviate from its linear path T2 = T1 + √2Ra in the region closer to the origin. It is recognized that two causes preclude (0) the state T2 = T1 being realized; on one hand the finite value of Ra, and on the other hand the non-vanishing value of δ that governs the orientation of the ”jump”. For larger values of δ, viz. a higher relative value of V1 , the slope is more accentuated and the slow-motion part becomes restricted to a neighborhood of the junction point C where the dynamics of the system are dominated by non-linear effects. For an analogous examination of the slow motion in the sector T22 < T12 , when Ra >> 1, the same approach as in Section 6 is taken. Hereby the zeroth-order problem assumes the form (0) 2

−(T22 − T1

(0)

)(T2 − T1 ) + (

θ (0) (0) − T1 ) + κ(T2 − T1 ) = 0, 1−θ (0)

where = 1/Ra. After expanding T1 (0)

T1

(0,0)

= T1

as

+ R(T2 ) + 2 S(T2 ) + ... , 33

(148)

(149)



+ ... (147)

insertion in the equation above and ordering in powers of yields (0,0)2

(0,0)

0 : (T22 − T1

1

3/2

)(T2 − T1 ) = 0, =⇒ 2κ + 1 θ/(1 − θ) : R(T2 ) = + , 4T2 4T22 R(T2 )2 (1 + κ)R(T2 ) : S(T2 ) = − . 4T2 4T22

(0,0)

T1

= −T2 ,

(150) (151) (152)

The structure of the higher-order results above implies that for increasing (0) values of T2 when Ra >> 1, T1 → −T2 in agreement with the limit cycle representations shown in Figs. 8 and 10. Note furthermore that here the firstorder correction term R(T2 ) incorporates the effects of both the parameters θ and κ. A further important property of the asymptotic analysis, which becomes particularly striking for Ra >> 1, is revealed when the rapid phases of motion dealt with in Sections 4 and 7 are examined. As an example, Eq. (120), dealing with the ”jump” initiated at the modified cusp point A∗ , can be reformulated as  T1  T1 (k4 − k2 )z − (k3 + k5 ) B T2 (T1 , δ) = T2 − δ dz − δ dz + ...(153) ∗ ∗ z 3 − k1 z 2 − k4 z + k5 T1A T1A Insights concerning this relationship can be gained by adding Eqs. (1) and (2): d 1+θ (T2 + δT1 ) = − T1 − T 2 , dt 1−θ

(154)

where T2 + δT1 , represents the total heat content of the system. Integrating this equation over the time τ required for the solution to move from the modified cusp point A∗ to the auxiliary drop point B † and comparing with † Eq. (153) evaluated at T1 = T1B yields the identity 

τ 0

1+θ − T1 − T2 )dt = −δ ( 1−θ



T1B

†

∗

T1A

(k4 − k2 )z − (k3 + k5 ) dz... . (155) z 3 − k1 z 2 − k4 z + k5

When the system is subjected to strong forcing, i.e. Ra >> 1, k4 → k2 , k5 → −k3 (cf. Section 7), Eq. (153) reduces to †

T2 (T1B , δ) + δT1B

†

∗

A∗

∼ T2A + δT1 ∼ 34

θ(δ − 1) . (1 + 2κ)(1 − θ)

(156)

For Ra >> 1 the heat content of the system is thus conserved in its entirety during this rapid phase of motion, i.e. the transition from A∗ to B † is adiabatic. This result is valid whatever the value of δ and only depends on the system being strongly forced and corresponds precisely to what previously has been shown by the present authors [21]. It is furthermore recognized that both of the integrals in Eq. (155) represent the deviation from a purely adiabatic transition. Consonant with our physical understanding of the system, Eq. (155) shows that strong forcing leads to a negligibly small duration τ of the ”jump”. It should be also underlined that for the limiting case δ → 0, and irrespective of the magnitude of the Rayleigh-type number Ra, the rapid phases of the discontinuous lowest-order solution takes place under adiabatic conditions. Note also that although the discussion above has limited itself to the ”jump” initiated at the modified cusp point A∗ , it is equally valid for the other rapid-motion segment of the limit cycle. For Ra >> 1, the asymptotic analysis to first order thus accommodates what in the theory of relaxation oscillations is known as a Mandelstam condition, cf. [26]. A constraint of this type, involving additional information, can be imposed on a lowest-order discontinuous solution. It is applied by prescribing that a system characteristic remains invariant during the rapid phases of a relaxation oscillation, in the present case T2 + δT1 , the total heat content of the system. These conditions are commonly based on independent considerations, but, as demonstrated above, the presently employed asymptotic formalism incorporates a Mandelstam-type condition to a lesser or greater degree, dependent upon the magnitude of Ra. This feature is in clear evidence from all of the phase-plane diagrams in this study, where the rapid-motion segments of the limit cycle manifest a slope determined by the magnitude of δ. The discussion in this section of solution properties for Ra >> 1 is of a certain general interest, since in these cases of strong forcing some of the physical mechanisms affecting the behavior of the system stand out particularly clearly and can be investigated analytically. The latter course of action is less feasible for smaller values of Ra; the strongly forced case is thus useful for illuminating processes which also are in effect for small values of Ra. 12. Alternative considerations Certain qualitative features of the system behavior are in particularly clear evidence at specific locations on the limit cycle. At the cusp near A, 35

where T1 = −T2 , the advective exchange between the vessels ceases and the heat conduction through the external walls of the vessels is minimal, whereas the conductive heat flux through the common wall of the system has a maximum. The latter heat flux, however, assumes its minimum value at the “turning-point” adjacent to B, since this is where the limit cycle is in its closest proximity to the symmetric diagonal T1 = T2 . This turning-point, where dT1 /dt = 0, is furthermore characterized by a momentary balance between the advective and conductive heat fluxes affecting the smaller vessel, which concurrently experiences a maximum heat loss through its “outer” wall. (Note that for Ra >> 1, application of the Mandelstam condition yields that the turning-point locus can be determined explicitly from a cubic algebraic equation in T1 .) At the other turning point, where dT2 /dt = 0 in the vicinity of C, a momentary balance is in force between the advective and conductive heat fluxes acting on the fluid in the larger vessel, which here also experiences a maximum heat gain through its outer wall. As the solution crosses the antisymmetric diagonal during its subsequent rapid motion in the direction of the drop point D, the advective exchange between the vessels briefly ceases, but resurges and has a local maximum before the slow motion in the sector T22 < T12 commences, cf. Fig. 7. From this figure it is also recognized that the advective heat exchange assumes its global maximum during the rapid “jump” from the cusp. In the limiting case of Ra >> 1, application of the Mandelstam condition permits an explicit determination of where this overall extremum is realized, viz. at ∗

TA T1 = 1 (1 + 2δ); 3

∗

TA T2 = − 1 (3 − 2δ), 3

(157)

a location which is shifted towards the cusp for increasing values of δ. Another approach facilitating an understanding of the system is based on Eq. (154), which is useful when examining the relaxation oscillation from a thermodynamic standpoint.The sign of the right-hand member of this equation (below denoted Ψ) specifies the nature of the overall thermal exchange between the global system and its surroundings. In a thermodynamic sense, the oscillation can be decomposed into two distinct stages: for Ψ < 0 an exothermal phase when the global system releases heat to its surroundings and for Ψ > 0 an endothermal phase of considerably longer duration when heat is absorbed. Given that 1+θ − T1 − T2 << 1, the transitions between 1−θ these phases take place near the antisymmetrical diagonal T1 = −T2 . The 36

slow motion in the sector T22 < T12 associated with a gradually decreasing heat flux into the system (cf. Fig. 7) thus represents the main part of the endothermic phase of the oscillation when the heat absorbed by the larger vessel exceeds that released the smaller one. The magnitudes of the overall heat fluxes to/from the larger/smaller vessel continue to decrease until an instability is triggered at the cusp, cf. Section 10. During the subsequent rapid motion (when only a minor exchange of heat between the global system and its surroundings take place) the exothermal phase commences. The slow motion in the sector T22 > T12 which follows is stable and associated with a comparatively rapid decrease of the heat released by the system, cf. Fig. 7. This exothermal phase is initially dominated by the conductive heat loss from the smaller vessel. As this gradually diminishes, T2 ultimately assumes its minimum value and a loss of stability ensues, leading to the global system reverting to its endothermal phase via rapid motion in the direction of the drop point D. From a teleological standpoint it might thus be argued that the oscillatory behavior of the system reflects the interplay between its two contradictory ”urges”. One, linked to the almost antisymmetric external forcing, is to maximize the temperature difference between the two vessels, i.e. to minimize | T1ex − T1 | and | T2ex − T2 |. This serves as a background to the endothermal phase-point motion in the proximity of T1 = −T2 dominated by the conductive heat fluxes. The other ”aspiration” of the system is to equalize the vessel temperatures, viz. to realize a state T1 = T2 . As, however, discussed in the previous section, this is precluded by the non-vanishing value of δ and by Ra being finite in magnitude. During this exothermal phase in the sector T22 > T12 the dynamics of the system are controlled by the interaction between the conductive and advective heat fluxes, which leads to decreasing values of T1 and T2 followed by a loss of stability and the subsequent transition of the system to its endothermal phase. 13. Discussion and conclusion Amending the lowest-order analysis [17] of the relaxation oscillation governed by Eqs. (1) and (2) with asymptotically valid corrections led, as seen, to significantly improved phase-plane representations of the limit cycle. The discontinuous approximation, valid in the limit δ → 0, has a number of shortcomings due to it being based on a degenerate version of Eq. (1) which assumes a perfect balance between the advective and conductive heat 37

fluxes affecting the fluid in the smaller vessel. This simplification had the consequence that the lowest-order analysis could not accommodate some of the more subtle time-dependent aspects of solution behavior. One ostensible inconsistency associated with the lowest-order analysis is that the rapid-motion phases are characterized by T2 = const., which formally implies that F (T1 , T2 ) = G(T1 , T2 ) = 0 during the “jumps”. As a consequence T2 is in fact continuously dependent upon T1 since T2 = (1 + θ)/(1 − θ) − T1 . The corrected solutions derived in Section 4 and 7 were capable of resolving this contradictory feature. As discussed in Section 11, these solutions can be seen as constituted by two parts, one which corresponds to the Mandelstam condition that the total heat content of the system be conserved during the rapid phases, cf. the previously reported results [21]. The other, which becomes negligible for Ra >> 1, accounts for the deviations from this purely adiabatic “ideal” case. Other disconcerting features associated with the discontinuous phases of the lowest-order solution are the abrupt transitions at the drop points B and D and at the junction point C. Proceeding from the auxiliary quantities B † and D† , the asymptotically valid corrections to the lowest-order results regularized these problems at the drop points, cf. Sections 5 and 8. The comparatively subtle analysis carried out in Section 3 served the same purpose at the junction point C. As evident from the phase-plane representations of the relaxation oscillation in Figs. 4, 5, 8, and 10, the asymptotically valid corrections made these mergers take place in a smooth fashion. As regards the slow-motion segments of the limit cycle it should be noted that in the sector T2 2 < T1 2 , the discontinuous approximation F (T1 , T2 ) = 0 does near-perfect justice to the numerical solution even for small values of Ra; in principle, as discussed in Section 10, no correction term is required here. In the sector T2 2 > T1 2 this correspondence was, however, less pronounced, and, as discussed in in Section 9, for weak thermal forcing the option arose of pursuing the analysis up to second order. As seen in Figs. 8 and 10, for Ra >> 1 the discontinuous solution F (T1 , T2 ) = 0 yielded perfect fits to the numerical results in both sectors of the phase-plane. From a comprehensive point of view the usefulness of the asymptotically valid correction terms derives from the fact that they take into account the actual imbalance between the advective and conductive fluxes affecting the smaller vessel (instead of postulating a perfect balance, as was the case in the discontinuous analysis). These corrections are also explicitly dependent on Eq. (2) governing the evolution in time of T2 , a process which was neglected 38

in the discontinuous approach. The results of the study may be summarized by underlining that the usefulness of the asymptotically valid corrections when representing the relaxation oscillations of the thermosyphon has been demonstrated, as has their linkage to the physics underlying the oscillatory behavior of the system. Finally an overview has been provided of some specific inconsistencies associated with the discontinuous lowest-order solution and how these were ameliorated by the asymptotically valid corrections. Acknowledgements The authors thank Drs. Marcus L¨ofverstr¨om and Saeed Falahat for valuable assistance in the preparation of this manuscript. Fariba Bahrami wishes to express her gratitude to the International Meteorological Institute at Stockholm University for support facilitating the completion of this study.

39

Appendix A: Lowest-order junction and drop points The junction point A, i.e. the cusp where the two branches of F (T1 , T2 ) = 0 coalesce, is given by T1A = −T2A = θ/[(1 + 2κ)(1 − θ)], cf. [17]. B, the associated drop point (T1B , T2A ), is found by determining T1B from a third2 2 degree algebraic problem originating from the constraint F (T2A > T1B ) = 0:   2 θ θ κ + 1 2 3 (T B ) − T1B + + (T1B ) + (1 − θ)(1 + 2κ) 1 (1 − θ)(1 + 2κ) Ra  3 θ (1 + κ)θ + − = 0. (A.1) Ra(1 − θ)(1 + 2κ) (1 − θ)(1 + 2κ) The discriminant of this equation is negative definite, corresponding to the irreducible case. One real root T1B = −T1A is already known, and hence the problem reduces to solving a second-degree algebraic problem in T1B . The pertinent solution is T1B

θ − =− (1 − θ)(1 + 2κ)



κ+1 Ra

1/2 ,

(A.2) 2

2

the remaining root being spurious since we require that T1B > T2A . The junction point C is found from the condition ∂F (T1 , T2 )/∂T1 |T1C ,T2C = 0, yielding T1C

1 TC = 2 − 3 3



2 4(T2C )

κ+1 +3 Ra

1/2 ,

(A.3) 2

where the other root of the quadratic equation is disregarded since T2C > 2 T1C . The constraint F (T22 > T12 ) = 0 hereafter yields a quartic equation in T2C : θ 4κ2 + 44κ + 13 C 2 18(2κ − 1)θ C 3 (T2C ) + (T2 ) − (T ) + 1−θ 32Ra 32(1 − θ)Ra 2  2 θ (κ + 1)3 27 +4 − = 0. (A.4) 32(Ra)2 32Ra 1 − θ 4

(T2C ) −

Once T2C has been found, the drop point D at (T1D , T2C ) on the other branch of F (T1 , T2 ) = 0 is determined from the following third-degree algebraic 40

problem in T1D : 3 (T1D )

−

2 T2C (T1D )

 +

   θ κ+1 1 C 2 D C 3 C − (T2 ) T1 + (T2 ) − + κT2 = 0. (A.5) Ra Ra 1 − θ

Using equation (158), the discriminant Δ of this cubic equation becomes:  2 θ 64 C 4 64 θ 54 C 3 Δ= (T ) − (T ) + . (A.6) Ra 2 Ra 1 − θ 2 (Ra)2 1 − θ The branch F (T2 2 > T1 2 ) = 0 remains confined to the upper half-plane, and hence T2C > 0. Since furthermore −1 < θ < 0 (cf. Lundberg and Rahm [27]), Δ is positive definite, and thus the drop-point D is uniquely determined due to the third-degree algebraic problem in T2D only having one real root. Observe that for notational convenience, the representation of the lowestorder quantities dealt with in this appendix did not include the superscript (0). Appendix B: Equivalent integrals In Eq. (86) the integration is carried out in both sectors of the phaseplane, which in principle necessitates the use of the two functional forms of the integrand: H1 (x, T2C ) =

−L1 + L2 ; L 1 + L3

x > T1 ,

(B.1)

H2 (x, T2C ) =

L1 + L2 ; −L1 + L3

x < T1 ,

(B.2)

where T1 is determined from Eq. (88), viz T1 = −T2 (T1 , δ), and 2

L1 (x, T2C ) = Ra(T2C − x2 )(T2C − x), L2 (x, T2C ) =

1 − T2C − κ(T2C − x), 1−θ

L3 (x, T2C ) =

θ − x + κ(T2C − x). 1−θ 41

(B.3) (B.4)

(B.5)

The integral in Eq. (86), here denoted Φ, thus assumes the form 

T1

Φ= b2

−L1 + L2 dx + L1 + L3



T1 T1

L 1 + L2 dx, −L1 + L3

(B.6)

whereafter a rearrangement shows that these integrals can be recast as  T1  T1 L2 − L 1 L 2 + L3 dx − 2 L1 2 dx = I1 − 2I2 . (B.7) Φ = L1 − L23 b2 L 1 + L 3 T1 As noted in section 4, Φ grows beyond bounds as T1 → T1D and hence the integration is only taken to extend to T1 = b3 . For the parameter range under consideration and the values of b2 and b3 used when representing the approximated limit cycles, | I2 /I1 | proved to be on the order of 10−2 . I2 can thus be neglected when evaluating Eq. (86), which simplifies the task of (1) determining the first-order correction T2 (T1 ) for this rapid-motion segment of the limit cycle. Appendix C: Some useful calculations Note that for computational purposes all following expressions are to be evaluated at the junction point C. 





M Gξ FT2 − GFT2 ξ γ˜ξ = − ,  γ(0, 0) (FT2 )2 



M Gη γ˜η = −  , γ(0, 0) FT2 



 η2

γ˜

M Gη 2 =−  , γ(0, 0) FT2 



(3)



F T2 ξ  M Gξ2 FT2 − GFT2 ξ2 γ˜ = − − 2 γ˜ ,   γ(0, 0) (FT2 )2 F T2 ξ  ξ2

42





γ˜ηξ

(3)



γ˜ηξ2

γ˜

γ˜







(3)











(3)

(4)

M Gξ3 FT2 + Gξ2 FT2 ξ − Gξ FT2 ξ2 − GFT2 ξ3 =− +  γ(0, 0) (FT2 )2   2 (3)  F T2 ξ 2  F T2 ξ FT2 ξ   γξ , −2  γ˜ξ − 4  γ˜ξ2 − 2˜  F T2 F T2 F T2 (4)

(4) ξ4



FT2 ξ  M Gηξ2 FT2 − Gη FT2 ξ2 =− − 2 γ˜ ,   γ(0, 0) (FT2 )2 FT2 ηξ (3)

(3) ξ3



M Gηξ FT2 − Gη FT2 ξ =− ,  γ(0, 0) (FT2 )2

(3)





(4)

(5)

M Gξ4 + 2Gξ3 FT2 ξ − 2Gξ FT2 ξ3 − GFT2 ξ4 =− +  γ(0, 0) (FT2 )2   2 (3)  FT2 ξ  FT2 ξ 2  (3) FT2 ξ γξ2  − 6˜ γξ 2 −6˜ γξ3  − 6˜  F T2 F T2 F T2 



−6˜ γξ

(3)

F T2 ξ F T2 ξ 2 

FT22

(4)



− 2˜ γξ

F T2 ξ 3 

F T2

,

    Gξ (ξ, η) = M Ra −(T2C − η)2 − 2(T1C − M ξ)(T2C − η) + 3(T1C − M ξ)2 − κ , 

Gξ2 (ξ, η) = 2Ra.M 2 ((T2C − η) − 3(T1C − M ξ)), (3)

Gξ3 (ξ, η) = 6RaM 3 , (4)

Gξ4 (ξ, η) = 0, 

Gη (ξ, η) = −Ra((T1C − M ξ)2 + 2(T1C − M ξ)(T2C − η) − 3(T2C − η)2 ) + (1 + κ), 43



Gη2 (ξ, η) = 2Ra((T1C − M ξ) − 3(T2C − η)), 

Gηξ (ξ, η) = 2M Ra((T1C − M ξ) + (T2C − η)), (3)

Gηξ2 (ξ, η) = −2M 2 Ra,   FT2 (ξ, −ξ 2 ) = Ra −(T1C − M ξ)2 + 3(T2C + ξ 2 )2 − 2(T1C − M ξ)(T2C + ξ 2 ) + κ,   FT2 ξ (ξ, −ξ 2 ) = 2Ra (T1C − M ξ) (M − 2ξ) + (T2C + ξ 2 )(M + 6ξ) ,  (3) FT2 ξ2 (ξ, −ξ 2 ) = 2Ra −M (M − 2ξ) − 2(T1C − M ξ) + 2ξ(M + 6ξ) + 6(T2C + ξ 2 ) , (4)

FT2 ξ3 (ξ, −ξ 2 ) = 12Ra(M + 6ξ), (5)

FT2 ξ4 (ξ, −ξ 2 ) = 72M Ra.

44

Bibliography [1] C. F. Bonilla, Nuclear Engineering, (McGraw-Hill, New York 1957). [2] Y. Zvirin, A review of natural circulation loops in pressurized water reactors and other systems, Nucl. Eng. and Design 67(1981), 203-225. [3] I. S. Kyung and S. Y. Lee, Periodic flow excursion in an open two-phase natural circulation loop, Nucl. Eng. and Design 162(1996), 223-232. [4] C. Wu, S. Wang and C. Pan, Chaotic oscillations in a low pressure twophase circulation loop under low power and and high inlet subcooling conditions. Nucl. Eng. and Design 162(1996), 233-244. [5] G. Morrison and J. Braun, System modeling and operation characteristics of thermosyphon solar water heaters, Solar Energy 34(1985), 389405. [6] R. Misra, Thermal stratification with thermosyphon effects in solar water heating systems, Energy Conv. Mgmt 33(1994), 193-203. [7] M. Bernier and B. Balaga, A 1-D/2-D model and experimental results for a closed-loop thermosyphon with vertical heat transfer sections, Int. J. Heat Mass Transf. 35(1992), 2969 -2982. [8] A. Wagner, Review of thermosyphon applications, Cold. Reg. Res. Eng. Lab. Report ERDC/CRREL TR-14-1, 2014 [9] Y. Lai, H. Guo and Y. Dong, Laboratory investigation on the cooling effect of the embankment with L-shaped thermosyphon and crushedrock revetment in permafrost regions. Cold Reg. Sci. and Tech. 58(2009), 143-150. [10] A. Franco and S. Filippeschi, Closed loop two-phase thermosyphon of small dimensions: a review of the experimental results, Micrograv Sci. Technol. 24(2012), 165-179. [11] C. Chang, S. Kuo, M. Ke and S. Chen, Two-phase closed- loop thermosyphon for electronic cooling, Exp. Heat Transf. 23(2010), 144-156. [12] P. Welander, On the oscillatory instability of a differentially heated fluid loop, J. Fluid Mech. 29(1967), 17-30. 45

[13] H. Creveling, J. de Paz, J. Baladi and R. Schoenhals, Stabilty characteristics of a single-phase free convection loop, J. Fluid Mech. 67(1975), 65-84. [14] M. Sen, E. Ramos and C. Trevino, The toroidal thermosyphon with known heat flux, Int. J. Heat Mass Transf. 28(1985), 219-233. [15] L. Rubenfeld and W. Siegmann, Nonlinear dynamic theory for a doublediffusive convection model, SIAM J. Appl. Math. 32(1977), 871-894. [16] B. van der Pol, On ”Relaxation-Oscillations”, Phil. Mag., 2(1926), 978992. [17] P. Lundberg, F. Bahrami, and M. Zarroug, A note on the asymptotic analysis of a thermal relaxation oscillator, Z. Angew. Math. u. Mech. 89(2009), No.12, 995-1001. [18] J. Grasman, Asymptotic Methods for Relaxation Oscillations and Applications, (Springer-Verlag, New York 1987). [19] F. Verhulst, Methods and Applications of Singular Perturbations, (Springer-Verlag, New York 2000). [20] E. F. Mishchenko and N. Kh. Rozov, Differential Equations with small Parameters and Relaxation Oscillations, (Plenum Press, New York 1980). [21] F. Bahrami, M. Zarroug, and P. Lundberg, Analysis of a themosyphon using a Mandelstam condition, Eng. Comp. Mech. 169(EM1)(2016), 2939. [22] J. Haag, Examples concrets d’´etude asymptotique d’oscillations de relaxation, Ann. Sci. Ecole Normale Super. 61(1944), 73-117. [23] D. Dorodnicyn, Asymptotic solution of a van der Pols equation, Priklad. Math. i Mekh., 11(1947), pp. 313-328. (In Russian.) Amer. Math. Soc, Transl., 88(1953). (English translation.) [24] G. N. Watson, A Treatise on the Theory of Bessel Functions, (Cambridge University Press 1958).

46

[25] J. Hadamard, Lectures on Cauchy’s Problem in Linear Partial Differential Equations, (Yale University Press, New Haven 1923). [26] N. Minorsky, Nonlinear oscillations, (D. Van Nostrand, Princeton, New Jersey 1962). [27] P. Lundberg and L. Rahm, A nonlinear convective system with oscillatory behaviour for certain parameter regimes, J. Fluid Mech. 139(1984), 237-260.

47

Figure 1: Schematic picture of the thermosyphon constituted by two well-mixed vessels connected by tubes through a thermally conducting wall. The presence of the adjoining thicker insulated wall indicates that differential thermal forcing can be applied to the system by immersing the vessels in separate constant-temperature baths.

48

Figure 2: Phase-plane diagram for for Ra = 25, θ = −49/50, and κ = 0.05 showing the fundamentals of the discontinuous lowest-order representation of the limit cycle. The heavy solid lines represent the two branches of F (T1 , T2 ) = 0 which merge at the cusp A. The stability in their immediate neighbourhood is indicated by the arrows signifying phasepoint motion. The solution becomes unstable at the cusp A and at the junction point C and ”jumps” to the drop points B and D.

49

Auxiliary functions

5

0

-5

-10 -4

-3

-2

-1

0

1

2

3

4

5

u

Figure 3: The auxiliary functions used for the analysis of the central part of the junction section at C. Legend: μν1 (u) dash-dotted line, μ2 ν2 (u) dashed line, μ3 ν3 (u) stippled, μ4 ν4 (u) weak solid line. The graph additionally includes the dividing solution ν0 (u), also known as the Haag special function, in the form of a heavy solid line. The results pertain to the standard case δ = 0.1, θ = −13/14, Ra = 10 and, κ = 0.05. Note that μ is proportional to δ 1/3 .

50

0.5

0.4

T2

0.3

0.2

0.1

-0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

T1

Figure 4: Limit cycles for δ = 0.05, θ = −13/14, Ra = 10, and, κ = 0.05. The asymptotically valid approximation is shown by the solid line, which coincides almost perfectly with its numerically counterpart shown by the dash-dotted line. The dashed line represents the discontinuous lowest-order results.

0.5

0.4

T2

0.3

0.2

0.1

-0.5

-0.4

-0.3

-0.2

-0.1

0

0.1

T1

Figure 5: As in Fig. 4, but for δ = 0.1. A discrepancy between the asymptotically valid representation of the limit cycle and its numerically determined counterpart can be seen at the merger between the rapid- and slow-motion segments in the sector T22 − T12 .

51

0.5 0.4 0.3

Temperatures

0.2 0.1 0 −0.1 −0.2 −0.3 −0.4 −0.5

2

3

4

5

6

7

8

9

10

11

Time

Figure 6: Numerically determined evolution in time of T1 (dash-dotted line) and T2 (dashed line) for δ = 0.1, θ = −13/14, Ra = 10, and κ = 0.05. The graph also includes the concurrent time evolution of the total heat content of the system T2 + δT1 (solid line).

0.6

Heat fluxes

0.4

0.2

0

−0.2

−0.4

2

3

4

5

6

7

8

9

10

11

Time

Figure 7: Numerically determined evolution in time of the heat fluxes affecting the system: the magnitude of the advective heat flux through the connecting tubes (solid line), the externally forced conductive heat fluxes through the outer walls of the smaller and larger vessel (dash-dotted and dashed line, respectively), and the magnitude of the heat flux through the conductive internal wall separating the two vessels (stippled line). These results pertain to the case of δ = 0.1, θ = −13/14, Ra = 10, and κ = 0.05.

52

0.5

0.4

T2

0.3

0.2

0.1

0 −0.5

−0.4

−0.3

−0.2

−0.1

0 T1

0.1

0.2

0.3

0.4

0.5

Figure 8: The case when δ = 0.1, θ = −13/14, Ra = 250, and κ = 0.05. The asymptotically valid representation of the limit cycle coincides almost perfectly with its numerically determined couterpart.

Correction term

10

1

0.1

0.01 0

0.05

0.1 T1 - T1C

0.15

0.2

Figure 9: The asymptotically valid corrections T2 (1) (T1 − T1C ) to the lowest-order slowmotion segment in the sector T22 > T12 as functions of the distance T1 −T1C to the singularity at the junction point C. Results are shown for Ra = 10 (dashed-dotted line), Ra = 25 (solid line), Ra = 50 (dotted line), Ra = 100 (dashed line), and Ra = 250 (solid and dotted line) and pertain to the case θ = −13/14 and κ = 0.05. Note the logarithmic scaling of the vertical axis.

53

0.5

0.4

T2

0.3

0.2

0.1

0 −0.5

−0.4

−0.3

−0.2

−0.1

0 T1

0.1

0.2

0.3

0.4

0.5

Figure 10: As in Fig. 8, but for δ = 0.2. A discrepancy between the asymptotically valid representation of the limit cycle and its numerically determined counterpart is discernible at the merger between the rapid- and slow-motion segments in the sector T22 − T12 .

54