Analysis of a thermosyphon using a Mandelstam condition

Bahrami, Fariba; Zarroug, Moundheur; Lundberg, Peter

doi:10.1680/jencm.15.00011

The characteristics of a thermally forced connected-vessel thermosyphon operating in an oscillatory mode have been determined using analytical techniques, the outcome of which is compared with results obtained by numerical integration of the governing equations. From a previous investigation it was known that adequate phase–plane representations of the limit cycles associated with oscillations could be obtained if the vessel-volume ratio was sufficiently small. This study aims at demonstrating how this constraint on the vessel volumes can be relaxed by prescribing a Mandelstam condition, that is, by postulating that the total heat content of the system remains conserved during the rapid phases of the oscillation. It was concluded that incorporating this Mandelstam condition in the analysis had the highly beneficial consequence that good analytical results could be obtained for much larger values of the vessel-volume ratio than those previously permitted.

1. Introduction

Thermosyphons have been of technical interest ever since the advent of the internal-combustion engine, when cooling systems utilised this principle for maintaining cylinder blocks at optimal working temperatures. In more recent times, thermosyphons have been used in nuclear-reactor engineering as well as for cooling computer processors. Interest in thermosyphons is, however, not limited to the engineering community; the underlying conceptual framework has – for example, also been applied to problems within the geological sciences (Torrance, 1979).

During the last few decades, interest in thermosyphons has seen a revival; this is due to their usefulness when exploiting alternative energy sources such as solar heating (e.g. Zvirin et al., 1977). A related example of the latter type of application comes from Scandinavia, where widespread use is made of heat pumps feeding on water at the lake bottom, which almost invariably is at temperatures of around 4°C – that is, the temperature at which fresh water assumes its maximum density.

In a previous study (Lundberg et al., 2009), the present authors examined a thermosyphon with water as the working fluid and operating over a temperature range centred at 4°C, a system described in detail in the next section. Under certain parametric conditions, the temperatures in the two interconnected reservoirs constituting the thermosyphon proved to manifest periodic behaviour comprising slow as well as rapid changes of the state of the system. In the limiting case of a very small ratio between the two vessel volumes, these oscillations could be dealt with analytically. For small values of the ratio, this asymptotically valid approach yielded phase–plane representations of the limit cycles corresponding to oscillatory behaviour of the system, which were in very good agreement with the results obtained by numerical integrating of the governing equations using a fourth-order Runge–Kutta scheme.

In this study, it will be demonstrated how this requirement of a very small vessel-volume ratio (so as to resolve the problem using analytical techniques) can be relaxed by the imposition of a Mandelstam condition (Minorsky, 1962) which assumes that the system has a property which remains unchanged during the rapid phases of the oscillation. For the thermosyphon under consideration here, the quantity that may be regarded as conserved is, as will be shown, the total heat content of the system. Using a number of examples it is shown how the extended analysis incorporating a Mandelstam condition for the rapid phases of the oscillation yields good phase–plane results even for comparatively large values of the vessel-volume ratio.

2. System description and mathematical model

The thermosyphon under consideration here consists of two well-mixed reservoirs of equal height H, but of different volumes V₁ and V₂ (cf. Figure 1, where the symbolic propellers indicate well-mixed states). The vessels are connected at their lower and upper ends by circular tubes of radius r and length l/2 through a thermally conducting common interface, here denoted as the ‘inner’ wall. The working fluid is fresh water of viscosity ν. For mathematical convenience, all temperatures will be given on a linearly transformed scale where 0° corresponds to 4°C – that is, the temperature at which water assumes its maximum density ρ₀. Since the density ρ displays a parabolic behaviour with regard to temperature around 4°C (Azouni, 1983), the equation of state with regard to the ‘shifted’ temperature T*∈ (− 4, + 4) thus assumes the form ρ = ρ₀(1 − αT*²), where α is the coefficient of thermal expansion. The temperatures in vessels 1 and 2 are T*₁ and T*₂, respectively. The system is subjected to differential thermal forcing by immersing vessels 1 and 2 in separate baths with the prescribed forcing temperatures T₁^ex and T₂^ex, respectively. These constant-temperature baths give rise to heat transfer to or from the working fluid through the thermally conducting ‘outer’ vessel walls of truncated cylindrical shape.

Figure 1.

The physical system consists of two well-mixed (cf. the symbolic propellers) vessels of equal height H, where the small vessel is of volume V1 and the large one of volume V2. At their upper and lower ends, these two reservoirs are connected by tubes of radius r and lengths l/2 through a thermally conducting ‘inner’ wall (shaded). The conducting ‘outer’ vessel walls of truncated cylindrical shape permit heat transfer to or from the external baths of constant forcing temperatures T1ex and T2ex. The acceleration of gravity g acts downwards Refer to the image caption for details.

View large Download slide

The physical system consists of two well-mixed (cf. the symbolic propellers) vessels of equal height H, where the small vessel is of volume V₁ and the large one of volume V₂. At their upper and lower ends, these two reservoirs are connected by tubes of radius r and lengths l/2 through a thermally conducting ‘inner’ wall (shaded). The conducting ‘outer’ vessel walls of truncated cylindrical shape permit heat transfer to or from the external baths of constant forcing temperatures T₁^ex and T₂^ex. The acceleration of gravity g acts downwards

A density difference between the fluids in the well-mixed reservoirs drives a flow q – governed by the law of Hagen–Poiseulle – through the lower tube, the upper acting as a passive conduit

q = \frac{r^{4} g π α H}{8 l ν} (T_{2}^{* 2} - T_{1}^{* 2})

1

where g is the acceleration due to gravity. Heat transfer also takes place between the vessels by thermal conduction through their common ‘inner’ wall (shaded in Figure 1).

The rate of change of temperature in each of the vessels is characterised by the time scale τ and is governed by the heat fluxes affecting the vessel in question – that is, the advective fluxes through the connecting tubes and the conductive fluxes through its ‘outer’ and ‘inner’ walls. These considerations give rise to the following rescaled equations governing the evolution of the non-dimensionalised vessel temperatures T₁ = T₁^*/ΔT and T₂ = T₂^*/ΔT (where ΔT = T₁^ex − T₂^ex) in the small and large vessels, respectively

\begin{aligned} δ \frac{d T_{1}}{d t} = & Ra |T_{2}^{2} - T_{1}^{2}| (T_{2} - T_{1}) + \frac{1}{Pe} (\frac{θ}{1 - θ} - T_{1}) \\ + \frac{κ}{Pe} (T_{2} - T_{1}) = F (T_{1}, T_{2}) \end{aligned}

2

\begin{aligned} \frac{d T_{2}}{d t} = & - Ra |T_{2}^{2} - T_{1}^{2}| (T_{2} - T_{1}) + \frac{1}{Pe} (\frac{1}{1 - θ} - T_{2}) \\ - \frac{κ}{Pe} (T_{2} - T_{1}) = G (T_{1}, T_{2}) \end{aligned}

3

The thermal response of the system is governed by five dimensionless parameters arising from the scaling applied when non-dimensionalising the problem. $Ra = π τ α g {H r}^{4} (Δ T)^{2} / (8 {l ν V}_{2})$ ⁠, related to the classical Rayleigh number, reflects the strength of the thermal forcing; κ is the ratio between the ‘inner’- and ‘outer’-wall heat-transfer coefficients; $Pe = V_{2} / (κ τ)$ is a Péclet-type number inversely proportional to the heat-transfer coefficient ratio κ and the characteristic time scale τ; δ = V₁/V₂ is the vessel-volume ratio; and $θ = T_{1}^{e x} / T_{2}^{e x} < 0$ is the quotient between the prescribed forcing temperatures affecting the small and large vessels. The definitions of $Ra$ and $Pe$ deviate somewhat from those in standard use for the Rayleigh and Péclet numbers, and thus $Ra$ and $Pe$ as introduced above will in what follows be denoted Rayleigh- and Péclet-type numbers. For a detailed discussion of Equations 2 and 3, see Lundberg and Rahm (1984) as well as Lundberg (1996).

In a previous study (Lundberg et al., 2009), it was shown how oscillations of the system in the limit of a very small vessel-volume ratio, namely for δ → 0, could be described reasonably accurately in the (T₁, T₂) phase–plane by neglecting the left-hand side of Equation 2 – that is, by disregarding the time derivative $d T_{1} / d t$ ⁠, and basing the analysis on the stationary slow manifold of the system given by the non-analytical phase–plane curve F(T₁, T₂) = 0. These results are summarised in Figure 2, showing that the limit cycles corresponding to periodic behaviour of the temperatures T₁ and T₂ to lowest order mainly are constituted by motion on the two branches of the manifold. These ‘slow’ segments of the limit cycle are connected by rapid jumps, which in the limit of δ → 0 take place in parallel to the T₁-axis. On $F (T_{2}^{2} < T_{1}^{2}) = 0$ ⁠, the ‘left-hand’ branch of the manifold, the solution loses stability at the junction point A where the two branches of the manifold merge, which, as shown on the basis of algebraic considerations, takes place at

T_{1}^{A} = - T_{2}^{A} = \frac{θ}{{(1 + 2 κ) (1 - θ)}}

4

Figure 2.

The slow manifold F(T1, T2) = 0 of the system for Ra=25, Pe=1, θ = − 13/14 and κ = 0·05. The two branches of the manifold coalesce at A, their stability indicated by the arrows showing phase-point motion in their immediate neighbourhood. The solution loses stability at the junction points A and C, but resumes its stable motion along the slow manifold at drop points B and D Refer to the image caption for details.

View large Download slide

The slow manifold F(T₁, T₂) = 0 of the system for $Ra =$ 25, $Pe =$ 1, θ = − 13/14 and κ = 0·05. The two branches of the manifold coalesce at A, their stability indicated by the arrows showing phase-point motion in their immediate neighbourhood. The solution loses stability at the junction points A and C, but resumes its stable motion along the slow manifold at drop points B and D

On the ‘right-hand’ branch of the slow manifold – that is, on $F (T_{2}^{2} > T_{1}^{2}) = 0$ ⁠, instability sets in at the junction point C where T₂ has a minimum. Based on an application of the implicit-function theorem, the coordinates of C can be specified according to

T_{1}^{C} = \frac{T_{2}^{C}}{3} - \frac{1}{3} {(4 (T_{2}^{C})^{2} + 3 \frac{κ + 1}{Ra Pe})}^{1 / 2}

5

where $T_{2}^{C}$ is determined numerically using a Newton–Raphson scheme applied to the following fourth-order algebraic equation

\begin{aligned} (T_{2}^{C})^{4} - \frac{θ}{1 - θ} (T_{2}^{C})^{3} + \frac{4 k^{2} + 44 k + 13}{32 Ra Pe} (T_{2}^{C})^{2} \\ - \frac{18 (2 κ - 1) θ}{32 (1 - θ) Ra Pe} + 4 \frac{(κ + 1)^{3}}{32 (Ra Pe)^{2}} \\ - \frac{27}{32 Ra Pe} {(\frac{θ}{1 - θ})}^{2} = 0 \end{aligned}

6

After stability has been lost at the junction points A and C, rapid motion parallel to the T₁-axis ensues, and subsequently the phase-point resumes its stable motion on the other slow-manifold branch at the drop points B and D, respectively.

It is now a straightforward matter to construct the limit cycle in the (T₁, T₂) plane by graphing T₂ as a function of T₁ from $T_{1}^{A}$ to $T_{1}^{D}$ after resolving the algebraic equation $F (T_{1}^{2} > T_{2}^{2}) = 0$ numerically (in a similar way as was done above for Equation 6) and hereafter proceeding in the same manner from $T_{1}^{C}$ to $T_{1}^{\tilde{B}}$ based on the numerical resolution of $F (T_{1}^{2} < T_{2}^{2}) = 0$ ⁠. Since $T_{2}^{B} = T_{2}^{A}$ and $T_{2}^{D} = T_{2}^{C}$ ⁠, the junction points A and C are hereafter connected with the drop points B and D, respectively, by straight lines parallel to the T₁-axis, which completes the construction of the limit cycle for the limit δ → 0.

Note, however, that a fundamental feature of this lowest-order analysis is that the effects of a non-zero value of δ have been completely disregarded.

3. Application of a Mandelstam condition

The consequences of a finite vessel-volume ratio δ can, however, to a certain extent be taken into account by applying a Mandelstam condition (Minorsky, 1962). This can be done when the system is characterised by a quantity that remains unchanged during the rapid phases of the oscillation; here the total heat content of the two vessels is taken to be conserved during these transitions, which thus takes place under adiabatic conditions.

A formal rationale for the application of a condition of this type to the problem presently under consideration is obtained by adding the two governing equations, which yields

\frac{d ({δ T}_{1} + T_{2})}{d t} = \frac{1}{Pe} (\frac{θ + 1}{θ - 1} - (T_{1} + T_{2}))

7

By definition, the Péclet-type number $Pe$ is inversely proportional to the characteristic time scale. During sufficiently rapid transitions between junction and drop points, the non-dimensional heat content δT₁ + T₂ of the system thus remains approximately conserved; this is due to the fact that during the sufficiently rapid phases of the oscillation there is insufficient time for any heat transfer to take place through the ‘outer’ walls of the vessels which are exposed to the baths of constant dimensional temperatures $T_{1}^{ex}$ and $T_{2}^{ex}$ that force the system.

The junction points A and C are already known, cf. Equations 4–6, and thus the modified drop points $\tilde{B}$ and $\tilde{D}$ can be determined using analytical techniques.

$\tilde{B}$ is found from the condition that $(T_{1}^{\tilde{B}}, T_{2}^{\tilde{B}})$ is to satisfy $F (T_{2}^{{\tilde{B}}^{2}} > T_{1}^{{\tilde{B}}^{2}}) = 0$ under the constraint that the total heat content of the system remains unchanged during the transition from A to B, which latter condition yields

T_{1}^{\tilde{B}} = \frac{1}{δ} (\frac{θ (δ - 1)}{(1 - θ) (1 + 2 κ)} - T_{2}^{\tilde{B}})

8

The condition that the drop point be located on the slow manifold results in the following cubic algebraic equation from which $T_{2}^{\tilde{B}}$ can be determined

\begin{aligned} (T_{2}^{\tilde{B}})^{3} + \frac{θ}{1 - θ} \frac{(3 + 2 δ - δ^{2})}{(1 + 2 κ) (δ + 1)^{2}} (T_{2}^{\tilde{B}})^{2} + \frac{1}{(δ + 1)^{2}} (\frac{{κ δ}^{3} + (κ + 1) δ^{2}}{Ra Pe (δ - 1)} - \frac{θ^{2}}{(1 - θ)^{2}} \frac{2 δ + δ^{2} - 3}{(1 + 2 κ)^{2}}) T_{2}^{\tilde{B}} \\ + \frac{θ}{(1 - θ)} \frac{1}{(δ + 1)^{2}} (\frac{θ^{2}}{(1 - θ)^{2}} \frac{(δ - 1)^{2}}{(1 + 2 κ)^{3}} + \frac{δ^{2}}{Ra Pe} (\frac{δ}{δ - 1} + \frac{κ + 1}{1 + 2 κ})) = 0 \end{aligned}

9

Since one root $T_{2}^{\tilde{B}} = T_{2}^{A}$ corresponding to the cusp at A (located on the unstable part of the $F (T_{2}^{{\tilde{B}}^{2}} > T_{1}^{{\tilde{B}}^{2}}) = 0$ manifold branch, cf. Figure 2) is already known, the problem reduces to a quadratic equation with the solution

T_{2}^{\tilde{B}} = \frac{θ}{1 - θ} \frac{1}{1 + 2 κ} \frac{δ - 1}{δ + 1} + {(\frac{δ^{2} (κ δ + κ + 1)}{(1 - δ) (δ + 1)^{2} Ra Pe})}^{1 / 2}

10

with the other root being spurious since $T_{2}^{{\tilde{B}}^{2}} > T_{1}^{{\tilde{B}}^{2}}$ is required. $T_{1}^{\tilde{B}}$ has already been specified in terms of $T_{2}^{\tilde{B}}$ ⁠, and hence the modified drop point $\tilde{B}$ is now fully determined.

The other modified drop point $\tilde{D}$ is calculated by requiring $(T_{1}^{\tilde{D}}, T_{2}^{\tilde{D}})$ to satisfy $F (T_{2}^{{\tilde{D}}^{2}} < T_{1}^{{\tilde{D}}^{2}}) = 0$ under the adiabatic constraint

T_{1}^{\tilde{D}} = \frac{1}{δ} ({δ T}_{1}^{C} + T_{2}^{C} - T_{2}^{\tilde{D}})

11

The resulting cubic equation in $T_{2}^{\tilde{D}}$ is

\begin{aligned} Ra Pe (1 - \frac{1}{δ^{2}}) (1 + \frac{1}{δ}) (T_{2}^{\tilde{D}})^{3} + Ra Pe (1 + \frac{1}{δ}) (1 - \frac{3}{δ}) \\ \times {(T_{1}^{C} + \frac{T_{2}^{C}}{δ})}^{2} (T_{2}^{\tilde{D}})^{2} - (κ + 1) (T_{1}^{C} + \frac{T_{2}^{C}}{δ}) \\ + \frac{1}{δ} (κ + 1) + κ) T_{2}^{\tilde{D}} + \frac{θ}{1 - θ} - Ra Pe {(T_{1}^{C} + \frac{T_{2}^{C}}{δ})}^{3} \\ + Ra Pe (1 + \frac{3}{δ}) (T_{1}^{C} + \frac{T_{2}^{C}}{δ}) = 0 \end{aligned}

12

As shown in the Appendix, this third-degree algebraic equation has a unique real solution $T_{2}^{\tilde{D}}$ ⁠, and thus the modified drop point $\tilde{D}$ associated with the adiabatic transition from the junction point C is also fully determined.

Once the junction points A and C as well as the modified drop points $\tilde{B}$ and $\tilde{D}$ are known, the corresponding limit cycle can be constructed in the same way as previously was done in the limiting case when δ → 0, as reported in Section 2.

In the next section, the authors apply this analytical scheme for representing the limit cycle in the (T₁, T₂) plane for a number of cases based on the parameter regime dealt with in the previous investigations of Equations 2 and 3 (Lundberg, 1996; Lundberg and Rahm, 1984; Lundberg et al., 2009). This choice of parameters serves to facilitate direct comparisons with previous results concerning Equations 2 and 3.

4. Limit-cycle results

Figures 3–7 show a number of limit cycles in the (T₁, T₂) phase–plane, which have been assembled using the technique outlined in the previous section. The constructed approximations incorporating the Mandelstam condition are shown by the weak solid lines, and the two branches of the slow manifold are included as dotted curves. For comparative purposes, Equations 2 and 3 have also been integrated numerically using a fourth-order Runge–Kutta scheme, and these limit-cycle results are also included in the figures as heavy solid lines. The cases shown are primarily for $P e =$ 1, θ = − 13/14 and κ = 0·05; the graphs illustrate the effects of varying the vessel-volume ratio δ and the Rayleigh-type number $Ra$ mirroring the strength of the thermal forcing. Figure 6, however, shows the consequences of increasing the heat-transfer coefficient ratio κ.

Figure 3.

View large Download slide

Approximated limit cycle constructed with the aid of a Mandelstam condition (thin solid line) and its numerically calculated counterpart (heavy solid line). The two branches of the slow manifold F(T₁, T₂) = 0 are dotted. The parameter values are $Ra =$ 9, $Pe =$ 1, δ = 0·20, θ = − 13/14 and κ = 0·05. The rapid segments of the limit cycle show a pronounced slope induced by the condition that the total heat content of the system remains unchanged during these transitions, which terminate at the modified drop points $\tilde{B}$ and $\tilde{D}$ on the right- and left-hand branches of F(T₁, T₂) = 0, respectively

Figure 4.

View large Download slide

As in Figure 3, but for $Ra =$ 250. Note the increased size of the limit cycle

Figure 5.

View large Download slide

As in Figure 3, but for δ = 0·001 and $Ra =$ 1·673 – that is, a very small vessel-volume ratio and weak forcing. Note the aberrant scaling of the coordinate axes compared with Figures 3, 4, 6 and 7

Figure 6.

View large Download slide

As in Figure 4, but here for a large value of the heat-transfer coefficient ratio, viz. κ = 0·70. Note the pronounced decrease in size of the limit cycles compared with those in Figure 4

Figure 7.

View large Download slide

As in Figure 3, but for $Ra =$ 50 and a very large vessel-volume ratio δ = 0·50. Note the discrepancy between the constructed approximation of the limit cycle and its numerically determined counterpart

A comparison between Figures 3 and 4, both for a moderately large vessel-volume ratio δ = 0·20 but for $Ra =$ 9 and 250, respectively, shows the effects of an increased thermal forcing of the system. A large value of $Ra$ not only increases the amplitude of the oscillation (viz. causes the limit cycle to expand) but also, as becomes evident when the two figures are compared, plays an essential role for ensuring a good correspondence between the approximated and numerically calculated limit cycles.

The effects of a small volume ratio are shown in Figure 5, where $Ra =$ 1·673 and δ = 0·001. As to be expected from the approximate analytical technique applied here, this minute value of the vessel-volume ratio δ ensures a very good fit between the approximated and numerically integrated limit cycles. The rapid segments of the limit cycle are nearly parallel to the T₁-axis, hereby almost entirely obviating the need to apply a Mandelstam condition. Note that the weak thermal forcing reflected in the diminutive value of $Ra$ implies an oscillation of small amplitude as recognised from the axes scales, which differ markedly from those used in Figures 3, 4, 6 and 7.

The approximated and numerically determined limit cycles shown in Figure 6 pertain to the case of $Ra =$ 250 and δ = 0·20, and show the effects of increasing the heat-transfer coefficient ratio to κ = 0·70. These results should be compared with those in Figure 4 for the same values of $Ra$ ⁠, δ and $Pe$ ⁠, and for κ = 0·05. It is recognised that a large value of κ causes the amplitude of the oscillation to decrease, as manifested by the small size of the limit cycle in Figure 6. This is primarily caused by the junction point A moving along the antisymmetric diagonal T₂ = − T₁ towards the origin of the phase–plane as κ increases, cf. Equation 4 yielding the cusp at A where the two branches of the manifold merge. For strong forcing, as signified by $Ra = 250$ ⁠, the location of the junction point C is far less affected by changes of κ, as seen by comparing Figures 4 and 6 as well as by noting that in these cases Equations 5 and 6 are comparatively insensitive to changes in κ.

Figure 7 shows the case of R = 50 and δ = 0·50. This large vessel-volume ratio not only manifests itself in the overall orientation of the rapid-transition segments of the limit cycles, but also in the analytically approximated limit cycle deviating considerably from the one calculated numerically, as expected for this comparatively large value of δ. This discrepancy, which becomes more pronounced as δ increases, is inevitable, given the comparatively straightforward approach (based on a heavily simplified version of the governing equations) taken when constructing the approximated limit cycle. The only analytical way to overcome this limitation of requiring a small δ would be to conduct a full-scale investigation analogous to the one undertaken by Dorodnicyn (1947) for the van der Pol equation. This task, however, lies beyond the scope of this study, the modest aim of which is to provide some results which are useful for practical purposes.

5. Oscillation-period results

The approximate procedure incorporating a Mandelstam condition can also be used to determine the limit-cycle periods. Equation 3 is reformulated as $d t = d T_{2} / G (T_{1}, T_{2})$ ⁠, which is integrated along the limit cycle to yield the oscillation period Ω. This line integral is dominated by the slow-manifold contributions, and hence it can be restated, for numerical purposes, in a convenient form

\begin{aligned} Ω = \int_{T_{1}^{C}}^{T_{1}^{\tilde{B}}} \frac{Pe (T_{2}^{2} + 2 T_{1} T_{2} - 3 T_{1}^{2}) + (κ + 1) / Ra}{((1 + θ) / (1 - θ) - T_{1} - T_{2}) (3 T_{2}^{2} - 2 T_{1} T_{2} - T_{1}^{2} + κ / (Ra Pe))} d T_{1} + \int_{T_{2}^{\tilde{D}}}^{T_{2}^{A}} \frac{Pe}{(1 + θ) / (1 - θ) - T_{1} - T_{2}} d T_{2} \end{aligned}

13

This expression can be evaluated numerically using Simpson's rule for the integrals. This has been done for two distinct cases, both with $Pe =$ 1 and θ = − 13/14, and for decreasing values of the vessel-volume δ = 0·20, 0·10 and 0·01. In the first case, the heat-transfer coefficient ratio was kept fixed at κ = 0·05 while $Ra$ is allowed to vary over its entire range as it is associated with oscillatory behaviour of the system. The results in the form of dashed lines are shown in Figure 8, which also comprise the periods obtained by integrating the governing equations numerically (solid lines). As foreseen, the overall accuracy of the results obtained using the approximate procedure is enhanced as δ decreases. The second case examined gave rise to similar results as given in Figure 9, which shows approximated (dashed lines) and numerically determined (solid lines) limit-cycle periods $Ω$ as functions of the heat-transfer coefficient ratio κ over its oscillatory range when the Rayleigh-type number was kept fixed at $Ra =$ 250.

Figure 8.

View large Download slide

Oscillation periods against the Rayleigh-type number $Ra$ for δ = 0·20, 0·10 and 0·001 when the heat-transfer coefficient ratio is kept fixed at κ = 0·05. The solid lines represent the numerically determined results obtained by integrating the governing equations; the dashed lines depict the periods obtained using the approximative procedure incorporating the Mandelstam condition. The other parameter values were $Pe =$ 1, θ = − 13/14

Figure 9.

View large Download slide

Oscillation periods against the heat-transfer coefficient ratio κ for δ = 0·20, 0·10 and 0·001 when the Rayleigh-type number is kept fixed at $Ra =$ 250. The solid lines represent the numerically determined results; the dashed lines show the periods obtained using the approximative procedure incorporating the Mandelstam condition. The remaining parameter values were prescribed as $Pe =$ 1 and θ = − 13/14

The periods are related to the overall size of the limit cycles, which conveniently may be defined as the distance between the junction points, namely $L = ((T_{1}^{A} - T_{1}^{C})^{2} + (T_{2}^{A} - T_{2}^{C})^{2})^{1 / 2}$ ⁠. Figure 10 shows two examples of how L varies as a parameter is changed. In one case (the solid curve), $Ra$ is taken to vary while the heat transport coefficient ratio is kept fixed at κ = 0·05; in the other case (the dashed curve), κ varies as the Rayleigh-type number is kept fixed at $Ra =$ 250. In both cases, the ‘standard’ values $Pe =$ 1 and θ = −13/14 were prescribed (note that the locations of the junction points A and C are independent of δ). A comparison between the results in Figure 10 and the periods shown in Figures 8 and 9 reveals overall the same type of parametric dependence upon changes of $Ra$ and κ – that is, L and the period increase/decrease as $Ra$ and κ, respectively, increase.

Figure 10.

View large Download slide

The diagram incorporates the results of two different evaluations aimed at determining the parametric dependence of the length-scale L characterising the limit cycle. The solid line represents L as a function of $^{10} \log Ra$ (lower horizontal scale) when the heat-transfer coefficient ratio is kept fixed at κ = 0·05. The dashed line shows L as a function of κ (upper horizontal scale) when the Rayleigh-type number is kept fixed at $Ra =$ 250. The Péclet-type number and the forcing-temperature ratio were prescribed as constant, viz. $Pe =$ 1 and θ = − 13/14. Note that the $Ra$ -dependent and the κ-dependent results show overall the same type of parametric dependence as manifested by the limit-cycle periods in Figures 8 and 9, respectively

The analytical approximation incorporating the Mandelstam condition permits a degree of δ-dependence of the approximated periods, shown in Figures 8 and 9, but the fact that their ranges as well as those of the corresponding numerically determined periods appear to be comparatively insensitive to changes of δ merits a comment. The reason is that in Figures 8 and 9 only rather small values of this vessel-volume ratio are taken into account. This implies that the directions of the ‘Mandelstam transitions’ from the junction point A to the modified drop point $\tilde{B}$ on the $F (T_{2}^{2} > T_{1}^{2}) = 0$ branch of the slow manifold only deviate moderately from the direction of the T₁-axis; consequently the residence time of the phase point on this branch is not much affected. As a contrast it may be instructive to examine Figure 7, pertaining to δ = 0·50, which shows that in this case the transition, orientated at a considerable angle from the T₁-axis, leads to the limit cycle by only encompassing a rather limited segment of the $F (T_{2}^{2} > T_{1}^{2}) = 0$ manifold branch. This has evident consequences, and a set of numerical calculations for the case of δ = 0·50, $Pe =$ 1, θ = − 13/14 and κ = 0·05 showed that the oscillation periods encompassed an interval from 3·51 (for $Ra$ = 12·09) to 5·03 (for $Ra$ = 91·6) – a considerable reduction compared with the period ranges visible in Figure 8. Note that due to the asymmetry of the slow manifold (Figure 2), the corresponding jumps from the junction point C to the modified drop point $\tilde{D}$ on the $F (T_{2}^{2} < T_{1}^{2}) = 0$ branch of the slow manifold are considerably shorter and hence associated with smaller changes of T₂. Consequently, the sensitivity to δ of the duration of the phase-point motion is less pronounced on this branch of the slow manifold.

A peculiarity of the numerically calculated oscillation periods, particularly prominent in Figure 9, is the period lengthening in immediate proximity of the parametric limits of $Ra$ and κ where the oscillatory solutions cease to exist. This is known as the ‘canard’ phenomenon (Lutz and Goze, 1981), and is caused by the phase point lingering on the unstable segment of the $F (T_{2}^{2} > T_{1}^{2}) = 0$ branch of the manifold before finally making the rapid transition to the other branch. This feature is intimately related to the fact that the parameter-dependent transitions of the thermosyphon from steady to oscillatory behaviour take place in the form of Hopf bifurcations (cf. Lundberg, 1996). Here rather subtle processes come into play, and the approximated periods reported in this study, that are based on a considerable simplification of Equations 2 and 3, cannot be expected to comprise this feature.

6. Discussion and conclusions

From an engineering perspective, the technologically dictated parameter regime ensuring a steady-state circulation of a thermosyphon is, for obvious reasons, of pre-eminent interest. For the system under consideration here, this question has been dealt with in previous studies (Lundberg, 1996; Lundberg and Rahm, 1984); hence, the present investigation has focused on dealing with the anomalous oscillatory behaviour that can occur under specific conditions and which may interfere with the smooth operation of the installation under consideration.

Given that simplified results can be useful in many situations, this study proceeded from the fact that it was known that the oscillatory behaviour of a thermosyphon could be reasonably well described using analytical techniques provided that the ratio of the vessel volumes was sufficiently small. This latter constraint, however, limited the practical application of the analytical procedure.

In this study, it has been demonstrated how to overcome this restriction by extending the previous analysis to also incorporate a Mandelstam condition based on the total heat content of the thermosyphon being conserved during the rapid phases of the oscillation. A series of comparisons, making use of the extended analytical procedure, was made between the analytically approximated and numerically determined phase–plane representations of the limit cycles associated with oscillatory behaviour of the system. This showed that a good correspondence between these two sets of results was obtained for even comparatively large values of the vessel-volume ratio. The same held true when the similarly determined periods of the oscillation were compared.

Thus, the conclusion of this study is that extending the previously established technique for analysing the thermosyphon (which only is valid in the limiting case of a very small vessel-volume ratio) by the application of a Mandelstam condition represents a significant improvement in that it permits the analysis to also be useful in cases of comparatively small volume differences between the reservoirs constituting the thermosyphon.

Appendix

A1.1 Unique real solution of the cubic equation

The discriminant of Equation 12 governing $T_{2}^{\tilde{D}}$ is

\begin{aligned} Δ = & \frac{δ^{6}}{108 (1 - δ^{2})^{6}} (27 \{δ^{2} (27 - 11 δ) {(T_{1}^{C} + \frac{T_{2}^{C}}{δ})}^{3} . + \frac{(1 - δ) δ}{3 Ra Pe} (T_{1}^{C} + \frac{T_{2}^{C}}{δ}) (κ δ + 1 - 2 κ) \\ {+ \frac{1 - δ}{Ra Pe} \frac{|θ|}{1 + |θ|}\}}^{2} + 4 {\{- \frac{4}{3} δ^{2} {(T_{1}^{C} + \frac{T_{2}^{C}}{δ})}^{2} + \frac{(κ + 1 + κ δ) (1 - δ)}{Ra Pe}\}}^{3}) \end{aligned}

14

where a prerequisite for oscillatory behaviour of the system is that θ be negative (Lundberg, 1996).

The expanded version of the cubic part of Δ comprises two negative terms

\begin{aligned} Q_{1} = & - \frac{256}{27} δ^{6} {(T_{1}^{C} + \frac{T_{2}^{C}}{δ})}^{6}; \\ Q_{2} = & - \frac{16 (κ + 1 + κ δ)^{2} (1 - δ)^{2} δ^{2}}{(Ra Pe)^{2}} {(T_{1}^{C} + \frac{T_{2}^{C}}{δ})}^{2} \end{aligned}

15

The quadratic part of the discriminant is only constituted by positive definite terms, and the ratio between |Q₁| and the square of the first term in the quadratic expression is

S_{1} = \frac{256}{729} \frac{δ^{2}}{27 - 11 δ}

16

which is unconditionally smaller than one since 0 < δ < 1. The ratio between |Q₂| and the square of the third term in the quadratic expression is

S_{2} = \frac{16}{27} \frac{(1 + |θ|)^{2}}{{|θ|}^{2}} (κ + 1 + κ δ)^{2} ({δ T}_{1}^{C} + T_{2}^{C})^{2}

17

The junction-point coordinates $(T_{1}^{C}, T_{2}^{C})$ are known from the analysis provided in Section 2. Using for the authors’ purposes a weak version of Equation 5, namely $T_{1}^{C} \approx - T_{2}^{C} / 3$ ⁠, it is possible to establish rough estimates of S₂. Given that θ = − 13/14 and that $T_{2}^{C}$ is postulated as assuming its definitional maximum value of 0·5, it is found that as δ ranges from 0·001 to the unrealistic value of 1, the maximum permissible κ-value compatible with S₂ < 1 increases from around 0·25 to 0·42 (the analogous results for θ = − 4/5 span the interval from κ ≈0·15 to 0·34, whereas for θ = −49/50, the maximum values range from κ ≈0·28 to 1·29). Under these highly conservative assumptions, S₂ is always smaller than one and since the same holds true for S₁, the discriminant Δ is definitely positive for the parameter ranges dealt with in this study. Hence, $(T_{1}^{\tilde{D}}, T_{2}^{\tilde{D}})$ is uniquely defined in these cases since the cubic equation has one real and two complex-conjugate roots.

Acknowledgements

The authors express their gratitude to not only two anonymous reviewers for highly constructive comments on a previous version of the manuscript which led to significant improvements, but also to Dr Saeed Falahat for valuable assistance in the preparation of this manuscript. Furthermore, Fariba Bahrami thanks the International Meteorological Institute at Stockholm University for support in facilitating the completion of this study.

Notation

$F (T_{1}, T_{2})$: right-hand side functional in Equation 1
$G (T_{1}, T_{2})$: right-hand side functional in Equation 2
g: acceleration due to gravity
$H$: vessel height
$L$: characteristic size of limit cycle
l/2: length of connecting tubes
$Pe$: Péclet-type number
$Q_{1}, Q_{2}, S_{1}, S_{2}$: auxiliary algebraic quantities used for estimating sign of $Δ$
q: volume flux based on law of Hagen–Poiseulle
$Ra$: Rayleigh-type number proportional to $(Δ T)^{2}$
r: radius of connecting tubes
T₁: non-dimensional water temperature in vessel 1
T₂: non-dimensional water temperature in vessel 2
$T_{1}^{*}$: dimensional water temperature in vessel 1
$T_{2}^{*}$: dimensional water temperature in vessel 2
(⁠ $T_{1}^{A}, T_{2}^{A}$ ⁠): coordinates of the junction point A
(⁠ $T_{1}^{B}, T_{2}^{B}$ ⁠): coordinates of the drop point B
(⁠ $T_{1}^{C}, T_{2}^{C}$ ⁠): coordinates of the junction point C
(⁠ $T_{1}^{D}, T_{2}^{D}$ ⁠): coordinates of the drop point D
(⁠ $T_{1}^{\tilde{B}}, T_{2}^{\tilde{B}}$ ⁠): coordinates of the modified drop point $\tilde{B}$
(⁠ $T_{1}^{\tilde{D}}, T_{2}^{\tilde{D}}$ ⁠): coordinates of the modified drop point $\tilde{D}$
$T_{1}^{ex}, T_{2}^{ex}$: temperatures of the external forcing bath
V₁, V₂: vessel volumes
α: thermal expansion coefficient of water
$Δ$: discriminant of cubic equation for $T_{2}^{\tilde{D}}$
$Δ T = T_{1}^{ex} - T_{2}^{ex}$: temperature scale
$δ = V_{1} / V_{2}$: vessel-volume ratio
$θ = T_{1}^{ex} / T_{2}^{ex}$: ratio between external forcing temperatures
κ: ratio between ‘inner’ and ‘outer’ wall heat-transfer coefficients
ν: water viscosity
ρ₀: maximum density of water
τ: characteristic time scale
Ω: period of limit cycle

REFERENCES

Azouni

M

(

1983

)

Hysteresis loop in water between 0° and 4°C

.

Geophysical and Astrophysical Fluid Dynamics

24

(

2

):

137

–

142

.

Google Scholar

Crossref

Dorodnicyn

D

(

1947

)

Asymptotic solutions of a van der Pol equation

.

Prikladnaya Matematika i Mekhanika

11

(

3

):

313

–

328

.

Google Scholar

Lundberg

P

(

1996

) Behaviour of a nonlinear convective oscillator as modified by sub- and super-critically unstable Hopf bifurcations.

In

Waves and Nonlinear Processes in Hydrodynamics

(

Grue

J

,

Gjevik

B

and

Weber

JE

(

eds

)).

Kluwer

,

Amsterdam, the Netherlands

, pp.

267

–

278

.

Google Scholar

Crossref

Lundberg

P

and

Rahm

L

(

1984

)

A nonlinear convective system with oscillatory behaviour for certain parameter regimes

.

Journal of Fluid Mechanics

139

:

237

–

260

.

Google Scholar

Crossref

Lundberg

P

,

Bahrami

F

and

Zarroug

M

(

2009

)

A note on the asymptotic analysis of a thermal relaxation oscillator

.

Zeitschrift für Angewandte Mathematik und Mechanik

89

(

12

):

995

–

1001

.

Google Scholar

Crossref

Lutz

R

and

Goze

M

(

1981

)

Non-Standard Analysis

.

Lecture Notes in Mathematics, vol. 881

.

Springer-Verlag

,

Berlin, Germany

.

Google Scholar

Crossref

Minorsky

N

(

1962

)

Nonlinear Oscillations

.

Van Nostrand

,

Princeton, NJ, USA

.

Google Scholar

Torrance

K

(

1979

)

Open-loop thermosyphons with geological applications

.

Journal of Heat Transfer

101

(

4

):

677

–

683

.

Google Scholar

Crossref

Zvirin

Y

,

Schitzer

A

and

Grossman

G

(

1977

)

The natural circulation solar heater-models with linear and nonlinear temperature distributions

.

International Journal of Heat and Mass Transfer

20

(

1

):

997

–

999

.

Google Scholar

Published with permission by the ICE under the CC-BY licence http://creativecommons.org/licenses/by/4.0/

2015

This is an open-access article distributed under the terms of the Creative Commons Attribution License, which permits unrestricted use, distribution, and reproduction in any medium, provided the original work is properly cited.

Analysis of a thermosyphon using a Mandelstam condition

1. Introduction

2. System description and mathematical model

3. Application of a Mandelstam condition

4. Limit-cycle results

5. Oscillation-period results

6. Discussion and conclusions

Appendix

A1.1 Unique real solution of the cubic equation

Acknowledgements

Notation

REFERENCES

Email Alerts

Cited By

Analysis of a thermosyphon using a Mandelstam condition Open Access

1. Introduction

2. System description and mathematical model

3. Application of a Mandelstam condition

4. Limit-cycle results

5. Oscillation-period results

6. Discussion and conclusions

Appendix

A1.1 Unique real solution of the cubic equation

Acknowledgements

Notation

REFERENCES

Email Alerts

Suggested Reading

Related Chapters

Recommended for you

Cited By

Analysis of a thermosyphon using a Mandelstam condition