Y B Fu

S P Pearce

K K Liu

Post-bifurcation analysis of a thin-walled hyperelastic tube under inflation

We consider the problem of bulging, or necking, of an infinite thin-walled hyperelastic tube that is inflated by an internal pressure, with the axial stretch at infinity maintained at unity. We present a simple procedure that can be used to derive the bifurcation condition and to determine the near-critical behaviour analytically. It is shown that there is a bifurcation with zero mode number and that the associated axial variation of near-critical bifurcated configurations is governed by a first-order differential equation that admits a locally bulging or necking solution. This result suggests that the corresponding bifurcation pressure can be identified with the so-called initiation pressure which featured in recent experimental studies. This is supported by good agreement between our theoretical predications and one set of experimental data. It is also shown that the Gent material model can support both bulging and necking solutions whereas the Varga and Ogden material models can only support bulging solutions. Relevance of the present method to the study of nonlinear wave propagation in a fluid-filled distensible tube is also discussed.

Introduction

Inflation of membrane tubes and the associated problem of bifurcation and instability is a classical subject that has been studied in many books and papers. For a membrane tube sealed by a rigid plug at each end and modeled by the Mooney-Rivlin material model, Kydoniefs and Spencer [START_REF] Kydoniefs | Finite axisymmetric deformations of an initially cylindrical elastic membrane[END_REF] showed that as the radius of the middle section increased, the pressure passed through a maximum (the limiting pressure) and then decreased, a phenomenon that had previously been observed for spherical membranes by Green and Shield [START_REF] Green | Finite elastic deformation of incompressible isotropic bodies[END_REF], Adkins and Rivlin [START_REF] Adkins | Large elastic deformations of isotropic materials. IX. The deformation of thin shells[END_REF], and Green and Adkins [START_REF] Green | Large elastic deformations and non-linear continuum mechanics[END_REF]. Kydoniefs and Spencer [START_REF] Kydoniefs | Finite axisymmetric deformations of an initially cylindrical elastic membrane[END_REF] also observed the tendency, in the case of a neo-Hookean membrane tube, for a nearly spherical bulge to form at the centre of the cylinder. Determination of the limiting pressure for other material models has also received some attention; see, for instance, Benedict et al [START_REF] Benedict | The determination of limiting pressure in simultaneous elongation and inflation of nonlinear elastic tubes[END_REF] and Alexander [START_REF] Alexander | The tensile instability of an inflated cylindrical membrane as affected by an axial load[END_REF]. Yin

A c c e p t e d m

a n u s c r i p t [START_REF] Yin | Non-uniform inflation of a cylindrical elastic membrane and direct determination of the strain energy function[END_REF] addressed the question whether a material could support a fully nonlinear bifurcation mode consisting of two uniform sections joined by a non-uniform transitional section, a phenomenon commonly observed in party balloons. He found the conditions that needed to be imposed on the strain-energy function for the material. The mechanisms for the initiation and propagation of bulges in an inflated membrane tube was made more transparent through the analysis of Chater and Hutchinson [START_REF] Chater | On the propagation of bulges and buckles[END_REF], and the experimental and numerical studies of Kyriakides and Chang [START_REF] Kyriakides | On the inflation of a long elastic tube in the presence of axial load[END_REF][START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF]. In particular, these authors observed that this problem shared the same features as a family of other problems such as propagating buckles in long metal tubes under external pressure [START_REF] Kyriakides | Experimental determination of the propagation pressure of circular pipes[END_REF][START_REF] Kyriakides | Propagating buckles in long confined cylindrical shells[END_REF], propagating necks in some polymeric materials when pulled in tension [START_REF] Hutchinson | Neck propagation[END_REF], and phase transition problems [START_REF] Ericksen | Equilibrium of bars[END_REF][START_REF] Bhattacharya | Microstructure of martensite: why it forms and how it gives rise to the shape-memory effect[END_REF][START_REF] Fu | Characterization and stability of two-phase piecewise-homogeneous deformations[END_REF]. To this family, we may also add the problem of kink band formation and propagation in layered structures [START_REF] Wadee | Kink band instability in layered structures[END_REF] and fibre-reinforced composites [START_REF] Fu | Continuum-mechanical modelling of kink-band formation in fibre-reinforced composites[END_REF]. We observe, however, that whereas the determination of the propagation pressure through the equal-area rule is now universally accepted, there is still some uncertainty concerning the determination of the pressure, termed the initiation pressure, at which bulging first occurs. Kyriakides and Chang [START_REF] Kyriakides | On the inflation of a long elastic tube in the presence of axial load[END_REF] took this pressure to be the limiting pressure corresponding to the first turning point in the pressure-volume diagram for uniform inflation, but this value was found to be slightly higher than what had been observed in the experiments (see their figure 9). Various other numerical schemes have also been used to determine the pressure-volume curves beyond the initiation pressure; see, for instance, Duffet and Reddy [START_REF] Duffett | The solution of multi-parameter systems of equations with application to problems in nonlinear elasticity[END_REF], Shi and Moita [START_REF] Shi | The post-critical analysis of axisymmetric hyper-elastic membranes by the finite element method[END_REF], Verron and Marckmann [START_REF] Verron | An axisymmetric B-spline model for the non-linear inflation of rubberlike membrane[END_REF], and the references therein. Additional experimental results are given recently by Pamplona et al [START_REF] Pamplona | Finite deformations of cylindrical membrane under internal pressure[END_REF] and Goncalves et al [START_REF] Goncalves | Finite deformations of an initially stressed cylindrical shell under internal pressure[END_REF].

The bifurcation and stability properties of the uniform inflation state are expected to be relevant to the determination of the initiation pressure. A thorough bifurcation study of the uniform inflation state was carried out by Haughton and Ogden [START_REF] Haughton | Bifurcation of inflated circular cylinders of elastic material under axial loading. I. Membrane theory for thin-walled tubes[END_REF], with some additional results given more recently by Haughton [START_REF] Haughton | Elastic Membranes, chapter 7 in Nonlinear Elasticity: Theory and Applications[END_REF]. Stability of the uniform inflation state was studied by Corneliussen and Shield [START_REF] Corneliussen | Finite deformation of elastic membranes with application to the stability of an inflated and extended tube[END_REF], Shield [START_REF] Shield | On the stability of finitely deformed elastic membranes; Part II: Stability of inflated cylindrical and spherical membranes[END_REF], and Chen [START_REF] Chen | Stability and bifurcation of finite deformations of elastic cylindrical membranes -part I. stability analysis[END_REF] using the energy method. However, although Shield's stability results were discussed by Kyriakides and Chang [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF], the results were not found to be relevant to the determination of the initiation pressure. This is probably because neither the bifurcation nor the stability analysis mentioned above is able to describe correctly the bifurcation mode having zero wave number, and it is in fact this mode that is relevant to the determination of the initiation pressure. This special mode is not sinusoidal or constant -its correct variation can only be found from a near-critical nonlinear analysis (see, for instance, Fu [START_REF] Fu | Perturbation methods and nonlinear stability analysis[END_REF]). It will be shown that this mode is localized so that it can also be supported by a finite tube. In fact, this localized mode is governed by a first-order differential equation that is an integral of the classical Korteweg-de Vries evolution equation.

The Gent material model [START_REF] Gent | A new constitutive relation for rubber[END_REF] has recently emerged as a popular and possibly realistic model for rubber and bio-materials. This is a so-called finite-chain model that imposes

A c c e p t e d m a n u s c r i p t

a maximum stretch that the material (or more precisely the polymeric molecular chains composing the material) can experience. This model has been discussed in a number of recent papers; see Horgan and Saccomandi [START_REF] Horgan | A description of arterial wall mechanics using limiting chain extensibility constitutive models[END_REF][START_REF] Horgan | Phenomenological hyperelastic strain-stiffening constitutive models for rubber[END_REF], and Ogden et al [START_REF] Ogden | On worm-like chain models within the threedimensional continuum mechanics framework[END_REF]. Gent [START_REF] Gent | Elastic instabilities of inflated rubber shells[END_REF] showed, for the case of closed ends, that the pressure-volume curve for the Gent material has a maximum followed by a minimum, which implies that for each pressure value in a certain interval there are two corresponding stable volume values (a property which is not found in the Varga or neo-Hookean material model). Gent [35, p.167] further stated that this is consistent with the experimental result that "thin-walled rubber tubes undergo a strikingly non-uniform deformation at a critical inflation pressure. One portion of the tube becomes highly distended as a bubble or aneurysm while the rest remains lightly inflated". This is exactly what has already been observed by Kyriakides and Chang [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF] for other material models.

The initial motivation for the current study is to extend Gent's [START_REF] Gent | Elastic instabilities in rubber[END_REF] own study and discover more about what the Gent model can predict with regard to the formation of an aneurysm (a localized, blood-filled dilation of a blood vessel caused by disease or weakening of the vessel wall). We focus on the simplest case when an infinite cylindrical membrane is inflated by an internal pressure with the remote axial stretch maintained at unity all the time (i.e. both before and after any bulging has taken place; this case differs from the closed-end case in that the pressure-volume curve for the Gent material would now be monotonically increasing). In the course of our study we find that analysis can be carried out much further than previously known before having to resort to numerical calculations.

The rest of this paper proceeds as follows. After formulating the problem in the next section, we show in Section 3 that the bifurcation condition can in fact be deduced from two simple equations that determine the two principal stretches at a point where the radius is a local maximum or minimum. In Section 4 we carry out a near-critical post-bifurcation analysis and derive the leading-order amplitude equation for the bifurcation solution. The amplitude equation is then used to determine under what conditions a localized solution is possible and, if it is, whether it is of the bulging or necking type. In Section 5 we present some numerical results for the fully nonlinear post-bifurcation states, performed with the aid of Mathematica, and compare them with the analytical near-critical results obtained in Section 4. In the final section, we discuss implications of our findings to the determination of the initiation pressure, and the relevance of our methodology to the study of nonlinear wave propagation in fluid-filled distensible tubes.

Governing equations

We consider the problem of inflation of an infinite cylindrical membrane tube that is incompressible, isotropic, and hyperelastic. The tube is assumed to have uniform thickness H and uniform inner radius R before inflation. To simplify analysis, we assume that the remote axial stretch is maintained at unity all the time, but our results can easily be modified for the case when the remote axial stretch takes any other value. We also assume that when the tube is inflated by an internal pressure, the inflated configuration maintains axial symmetry, but its radius can be constant or varying along the axial direction. We shall refer to the configuration with a constant radius the primary state and all other possible configurations bifurcated states. Thus, in general, the axisymmetric deformed configuration may be described by

r = r(Z), θ = Θ, z = z(Z), (2.1)
where Z and z are the axial coordinates of a representative material particle before and after inflation, respectively, and r is the inner radius after inflation; see Fig. 1. Since the deformation is axially symmetric, the principal directions of stretch coincide with the lines of latitude, the meridian and the normal to the deformed surface. Thus, the principal stretches are given by

λ 1 = r R , λ 2 = √ r ′2 + z ′2 , λ 3 = h H , (2.2)
where the indices (1, 2, 3) are used for the latitudinal, meridional and normal directions respectively and the primes indicate differentiation with respect to Z.

The principal Cauchy stresses σ 1 , σ 2 , σ 3 in the deformed configuration for an incompressible material are given by

σ i = λ i W i -p, i = 1, 2, 3 (no summation), (2.3)
where W = W (λ 1 , λ 2 , λ 3 ) is the strain-energy function, W i = ∂W/∂λ i , and p is the pressure associated with the constraint of incompressibility; see, for instance, Ogden [START_REF] Ogden | Non-linear elastic deformations[END_REF]. Utilising the incompressibility constraint λ 1 λ 2 λ 3 = 1 and the membrane assumption σ 3 = 0 we find

σ i = λ i Ŵi , i = 1, 2 (no summation), (2.4)

A c c e p t e d m a n u s c r i p t

where Ŵ (λ

1 , λ 2 ) = W (λ 1 , λ 2 , λ -1 1 λ -1
2 ) and Ŵ1 = ∂ Ŵ /∂λ 1 etc. From the fact that the resultant in the Z-direction at any cross section (balancing the tension in the membrane and the inflation pressure) must be a constant, we obtain

Ŵ2 z ′ λ 2 - 1 2 P λ 2 1 R = C 2 , (2.5)
where P is the internal pressure scaled by H and C 2 is a constant. By considering equilibrium in the normal direction of an infinitesimal volume element, we obtain

σ 2 (r ′′ z ′ -r ′ z ′′ ) λ 3 2 - z ′ σ 1 rλ 2 + P λ 3 = 0. (2.6)
With the use of (2.4) and the differentiated form of (2.5), we may integrate (2.6) once to obtain Ŵ -

λ 2 Ŵ2 = C 1 , (2.7)
where C 1 is another constant of integration. This integral was first noted by Pipkin [START_REF] Pipkin | Integration of an equation in membranes theory[END_REF]. Equations (2.5) and (2.7) correspond to Haughton's [START_REF] Haughton | Elastic Membranes, chapter 7 in Nonlinear Elasticity: Theory and Applications[END_REF] equations (3.29) and (3.27) respectively.

Without loss of generality, we take R = 1 (which is equivalent to scaling all the length variables by R). To find the constants C 1 and C 2 we assume that the cylinder is uniform far away from the origin with a constant radius r ∞ :

lim Z→±∞ r(Z) = r ∞ , lim Z→±∞ z(Z) = Z,
and hence

r ′ → 0, z ′ → 1, λ 1 → r ∞ , λ 2 → 1 as Z → ±∞.
Thus, from the outset we focus on localized bifurcation solutions by excluding sinusoidal-type solutions. On taking the limit Z → ∞ in (2.5), (2.6) and (2.7), we obtain

C 2 = Ŵ (∞) 2 - 1 2 P r 2 ∞ , P = Ŵ (∞) 1 r ∞ , C 1 = Ŵ (∞) - Ŵ (∞) 2 , (2.8)
where the superscript (∞) signifies evaluation at

λ 1 = r ∞ , λ 2 = 1.
We shall consider three types of materials: Varga, Ogden and Gent, for which the strainenergy functions are given, respectively, by

W = 2µ(λ 1 + λ 2 + λ 3 -3),
(2.9)

W = 3 r=1 µ r (λ αr 1 + λ αr 2 + λ αr 3 -3)/α r , (2.10)
W = - 1 2 µJ m ln(1 - J 1 J m ), J 1 = λ 2 1 + λ 2 2 + λ 2 3 -3, (2.11)

A c c e p t e d m a n u s c r i p t

where µ is the shear modulus for infinitesimal deformations, J m is a material constant and

α 1 = 1.3, α 2 = 5.0, α 3 = -2.0, µ 1 = 1.491, µ 2 = 0.003, µ 3 = -0.023
are material constants given by Ogden [START_REF] Ogden | Large deformation isotropic elasticity: on the correlation of theory and experiment for incompressible rubberlike solids[END_REF] (see also Ogden [36]). For rubbers, a typical value for J m is 97.2 which corresponds to a maximum stretch of 10 in uniaxial tension [START_REF] Gent | Elastic instabilities of inflated rubber shells[END_REF].

Horgan and Saccomandi [START_REF] Horgan | A description of arterial wall mechanics using limiting chain extensibility constitutive models[END_REF] give, from experimental data, values of J m between 0.422 and 3.93 for human arteries. In the limit J m → ∞, (2.11) reduces to

W = 1 2 µ(λ 2 1 + λ 2 2 + λ 2 3 -3), (2.12)
the energy function for the neo-Hookean material.

For the problem under consideration, the volume of the tube is given by v = πλ 2 1 per unit length in the z-direction and so turning points of the (P, v)-curve coincide with those of the (P, λ 1 )-curve. It is easy to show with the aid of (2.8) 2 that the (P, λ 1 )-curve has a maximum at λ 1 = r ∞ = 1.7321 and no minimum for the Varga material, has a maximum at r ∞ = 1.9582 and a minimum at r ∞ = 2.8380 for the Ogden material and is monotone increasing for the Gent material (for all J m > 0).

Bifurcation condition

With C 1 and C 2 known, equations (2.5) and (2.7) are two first-order coupled nonlinear differential equations for r(Z) and z(Z). For a bulging/necking solution that is symmetric about the origin Z = 0, we impose the conditions r ′ (0) = 0, z(0) = 0.

(3.1)

In order to integrate (2.5) and (2.7) from Z = 0 towards Z → ∞ to obtain a bulging/necking solution, we need to know r(0) as well as z(0). To find r(0), we evaluate (2.5) and (2.7) at Z = 0 and obtain respectively

Ŵ2 (r 0 , z ′ 0 ) - Ŵ (∞) 1 2r ∞ (r 2 0 -r 2 ∞ ) - Ŵ (∞) 2 = 0, (3.2) Ŵ (r 0 , z ′ 0 ) -z ′ 0 Ŵ2 (r 0 , z ′ 0 ) -Ŵ (∞) + Ŵ (∞) 2 = 0, (3.3)
where r 0 = r(0), z ′ 0 = z ′ (0) ≥ 0. These two equations can be solved to find r 0 (and z ′ 0 ) as a function of r ∞ . Obviously r 0 = r ∞ is always a solution which corresponds to the uniform solution r ≡ r ∞ of (2.5) and (2.7). A necessary condition for a non-uniform bifurcated solution to exist is that (3.2) and (3.3) have a solution other than r 0 = r ∞ , z ′ 0 = 1. In Figures 2(a of r ∞ . As expected, such solutions only exist for certain values of r ∞ . Of course, existence of such a non-trivial value of r 0 does not guarantee that it corresponds to a localized solution of r(Z). This will be discussed in more details in Sections 4 and 5. Now suppose that after the trivial solution r 0 = r ∞ has been factorized out, (3.2) and (3.3) have a solution given by r 0 = h(r ∞ ).

(3.4)

Then any bifurcation values are given by the roots of r ∞ = h(r ∞ ), i.e. the intersection of the trivial solution with the nontrivial solution. This suggests a method to derive the bifurcation condition as follows.

We assume that the solution of (3.2) for z ′ 0 in terms of r 0 is given by

z ′ 0 = k(r 0 ), with k(r ∞ ) = 1. (3.5)
We may alternatively choose (3.3) for this purpose. Of course, in general an explicit expression for k(r 0 ) does not exist, but we may differentiate (3.2) once and twice with respect to r 0 and then evaluating the resulting expressions at r 0 = r ∞ to obtain

Ŵ (∞) 1 - Ŵ (∞) 12 - Ŵ (∞) 22 k ′ (r ∞ ) = 0, (3.6) Ŵ (∞) 112 + 2k ′ (r ∞ ) Ŵ (∞) 122 + (k ′ (r ∞ )) 2 Ŵ (∞) 222 + k ′′ (r ∞ ) Ŵ (∞) 22 - Ŵ (∞) 1 r ∞ = 0. (3.7)
These two equations can be solved to obtain an explicit expression for k ′ (r ∞ ) and k ′′ (r ∞ ).

Next, with the use of (3.5), the left hand side of the second equation (3.3) is a function of r 0 only and we know that it having r ∞ as a non-trivial root defines the bifurcation condition. Differentiating the left hand side of (3.3) with respect to r 0 and then evaluating at r 0 = r ∞ followed by the use of (3.6), we find that the left hand side becomes identically zero, indicating that the trivial solution r 0 = r ∞ is in fact a double root. Thus, near r 0 = r ∞ the left hand side of (3.3), as a function of r 0 , expands like (r 0 -r ∞ ) 2 S(r ∞ ) + O((r 0 -r ∞ ) 3 ), and the bifurcation condition corresponds to S(r ∞ ) = 0. Therefore, the bifurcation condition can be derived by differentiating (3.3) twice with respect to r 0 and then evaluating at r 0 = r ∞ . We then obtain for the bifurcation condition,

r ∞ ( Ŵ (∞) 1 - Ŵ (∞) 12 ) 2 + Ŵ (∞) 22 ( Ŵ (∞) 1 -r ∞ Ŵ (∞) 11 ) = 0, (3.8)
which agrees with Haughton's [START_REF] Haughton | Elastic Membranes, chapter 7 in Nonlinear Elasticity: Theory and Applications[END_REF] equation (3.50) (note that the only plus sign in his equation (3.43) should be a minus). With the use of (2.8) 2 we may rewrite the above bifurcation condition as

r ∞ ( Ŵ (∞) 1 - Ŵ (∞) 12 ) 2 -r 2 ∞ Ŵ (∞) 22 dP dr ∞ = 0. ( 3.9)
Assuming that Ŵ (∞)

22

> 0 (see discussion in Section 6), we see from (3.9) that a bifurcation can never take place on the section of the (P, r ∞ )-curve where dP/dr ∞ < 0 (Haughton and Ogden [START_REF] Haughton | Bifurcation of inflated circular cylinders of elastic material under axial loading. I. Membrane theory for thin-walled tubes[END_REF]). A likely scenario is that as inflation proceeds, a bifurcation will take place before the limiting pressure is reached.

To illustrate the above idea, we consider the Varga material defined by (2.9). Equation (3.3) can be solved explicitly to express z ′ 0 in terms of r 0 . Equation (3.2) then becomes

(r 0 -r ∞ ) 2 r 0 - 2(1 + r 2 ∞ ) r 3 ∞ = 0, (3.10)
and so the bifurcation condition is given by r ∞ = 2(1 + r 2 ∞ )/r 3 ∞ , which has a unique positive real root r ∞ = 1 + √ 3 ≈ 1.6529, smaller than the value 1.7321 of r ∞ corresponding to the turning point of the (P, r ∞ )-curve.

Near-critical post-bifurcation analysis

In this section we show that the bifurcation condition (3.8) also emerges naturally in our near-critical post-bifurcation analysis. The latter analysis also yields the conditions under which a localized bifurcation solution is possible. We first rewrite (2.7) and (2.5) as

f (λ 1 , λ 2 ) ≡ Ŵ -λ 2 Ŵ2 -C 1 = 0, (4.1)
g(λ 1 , λ 2 ) ≡ λ 2 Ŵ2 (C 2 + P 2 λ 2 1 ) = z ′ , (4.2)
which define the two functions f and g. We write

λ 1 = r ∞ + y, (4.3)
where y is assumed to be small. Correspondingly, we expect λ 2 to have the expansion

λ 2 = 1 + d 1 y + 1 2 d 2 y 2 + O(y 3 ), (4.4)
where the coefficients can be obtained by substituting (4.3) and (4.4) into (4.1), expanding the latter as a power series in y, and then equating the coefficients of y and y 2 . This can best be carried out with the aid of a symbolic manipulation package such as Mathematica. Likewise, by substituting (4.3) and (4.4) into (4.2) and expanding in terms of y, we obtain

z ′ = 1 + g 1 y + 1 2 g 2 y 2 + O(y 3 ), (4.5)
where

g 1 = dg dy y=0 , g 2 = d 2 g dy 2 y=0
.

A c c e p t e d m a n u s c r i p t

On substituting (4.4) and (4.5) into λ 2 = √ r ′2 + z ′2 , we obtain d 1 = g 1 to leading order, and

y ′2 = ω(r ∞ )y 2 + O(y 3 ) (4.6)
to the next order, where ω(r

∞ ) = d 2 -g 2 . (4.7)
In terms of the various derivatives of Ŵ , we have

ω(r ∞ ) = -1 r ∞ Ŵ (∞) 2 Ŵ (∞) 22 r ∞ ( Ŵ (∞) 1 - Ŵ (∞) 12 ) 2 + Ŵ (∞) 22 ( Ŵ (∞) 1 -r ∞ Ŵ (∞) 11 ) . (4.8)
The behaviour of solution of the differential equation (4.6) changes when ω(r ∞ ) = 0. This gives the bifurcation condition consistent with the earlier expression (3.8).

For the Varga and Gent materials, ω(r ∞ ) is given by

ω(r ∞ ) = 2 + 2r 2 ∞ -r 4 ∞ 2 (r ∞ -1) r 2 ∞ , (4.9)
and

ω(r ∞ ) = J m r 2 ∞ (r 8 ∞ -6r 4 ∞ -4r 2 ∞ -3) -(r 2 ∞ -1) 2 (r 8 ∞ + 6r 6 ∞ + 8r 4 ∞ -2r 2 ∞ -1) r 2 ∞ (r 2 ∞ -1)(r 6 ∞ -(J m + 1)r 4 ∞ -(3J m + 1)r 2 ∞ + 1) , (4.10)
respectively. The corresponding expression for the Ogden material is too long and is not written out here. In Figures 2(b), 3(b) and 4(b) we have shown ω(r ∞ ) for the three materials. We observe that the Varga and Ogden materials each predict one bifurcation value, given by r ∞ = 1.6529 and r ∞ = 1.6955, respectively, but the Gent material predicts two bifurcation values: 1.74261 and 5.0850 when J m = 30, and 1.6907 and 9.6523 when J m = 97.2. We also find that as J m decreases the two bifurcation values move towards each other and coalesce at r ∞ = 2.2545 when J m = 11.0693. When J m is decreased further no bifurcation values exist. Thus, the type of bifurcation addressed in the present paper cannot occur in human arteries for which J m ranges between 0.422 and 3.93. The above expansion (4.6) can in fact be extended to higher orders with the aid of Mathematica. We have, in place of (4.6),

y ′2 = ω(r ∞ ) y 2 + γ(r ∞ ) y 3 + O(y 4 ), (4.11)
where the expression for the coefficient γ(r ∞ ) is easily obtained using Mathematica but is not written out here for brevity.

In the near-critical post-bifurcation analysis, we write

r ∞ = r cr + ǫr 1 , (4.12)
where r cr is a root of the bifurcation condition ω(r ∞ ) = 0, ǫ is a positive small parameter and r 1 is a constant. On substituting (4.12) into (4.11) and expanding ω around r cr , we obtain y ′2 = ω ′ (r cr )ǫr 1 y 2 + γ(r cr )y 3 , (4.13)

A c c e p t e d m a n u s c r i p t

where we have neglected the higher-order terms. We assume that y ′ does not change sign for 0 ≤ Z < ∞. Then when r 0 > r ∞ (a bulging solution), we have y ′ < 0 and y > 0; whereas when r 0 < r ∞ (a necking solution), we have y ′ > 0 and y < 0. Thus, for both cases we obtain from (4.13) y ′ = -y ω ′ (r cr )ǫr 1 + γ(r cr )y. (4.14)

The amplitude equation (4.14) indicates that y is of order ǫ and that its variation takes place on an O(1/ √ ǫ) lengthscale.

We now consider the solution of the following four prototypal equations:

y ′ 1 = -ay 1 1 + by 1 , y ′ 2 = -ay 2 1 -by 2 , (4.15
)

y ′ 3 = -ay 3 -1 + by 3 , y ′ 4 = -ay 4 -1 -by 4 , (4.16)
where y 1 , y 2 , y 3 , y 4 are functions of Z and a, b are positive constants given by

a = |ω ′ (r cr )ǫr 1 |, b = |γ(r cr )|/a 2 .
The relevance of these equations to (4.14) is obvious and we have

y =          y 1 if ω ′ (r cr )(r ∞ -r cr ) > 0, γ(r cr ) > 0, y 2 if ω ′ (r cr )(r ∞ -r cr ) > 0, γ(r cr ) < 0, y 3 if ω ′ (r cr )(r ∞ -r cr ) < 0, γ(r cr ) > 0, y 4 if ω ′ (r cr )(r ∞ -r cr ) < 0, γ(r cr ) < 0.
(4.17)

It can be shown that

y 1 = 1 b [sinh( 1 2 aZ + A)] -2 , if y 1 > 0 -1 b [cosh( 1 2 aZ + A)] -2 , if y 1 < 0 (4.18) y 2 = -1 b [sinh( 1 2 aZ + A)] -2 , if y 2 < 0 1 b [cosh( 1 2 aZ + A)] -2 , if y 2 > 0 (4.19) y 3 = 1 b[cos( 1 2 aZ + A)] 2 , y 4 = - 1 b[cos( 1 2 aZ + A)] 2 , (4.20)
where A is an integration constant which can be set to zero since the origin of Z can be arbitrarily chosen. By inspection, we may immediately dismiss y 3 , y 4 , y 1 > 0 and y 2 < 0 as we are only looking for localized solutions and they all blow up at finite values of Z. We are then left with

y = -1 b [cosh( 1 2 aZ)] -2 , if ω ′ (r cr )(r ∞ -r cr ) > 0, γ(r cr ) > 0, 1 b [cosh( 1 2 aZ)] -2 , if ω ′ (r cr )(r ∞ -r cr ) > 0, γ(r cr ) < 0. (4.21)
It follows that near-critical localized solutions exist only if

ω ′ (r cr )(r ∞ -r cr ) > 0, (4.22)

A c c e p t e d m a n u s c r i p t

and that when such a solution exists it is of the bulging type if γ(r cr ) < 0 and of the necking type if γ(r cr ) > 0. Thus, only a super-critical localized solution is possible if

ω ′ (r cr ) > 0, (4.23)
that is only if the curve of ω(r ∞ ) at the bifurcation point has a positive slope. Similarly, only a sub-critical localized solution is possible ω ′ (r cr ) < 0. When the condition (4.22) is not satisfied, our analysis above would predict a solution that blows up at a finite value of Z, but we should bear in mind the fact that our analysis has been conducted under the small-amplitude assumption so that terms of order y 4 in (4.11) could be neglected. Thus, when y grows and reaches a sufficient amplitude the higher order terms may no longer be neglected and they may act to suppress further growth. Our numerical integration to be explained in Section 5 does indeed show that whenever (4.22) is not satisfied the corresponding solution is usually periodic and hence does not satisfy our conditions at infinity.

We note that y and r 0 -r ∞ have the same sign and that a bulging/aneurysm solution corresponds to y > 0 whereas a necking solution corresponds to y < 0. Thus, referring to Figures 2(a), 3(a) and 4(a), we see that for the Gent material with J m = 97.2, near the first bifurcation point, r ∞ = 1.6907, the post-critical solution must necessarily be of the bulging type, whereas near the second bifurcation point, r ∞ = 9.6523, the post-critical solution must necessarily be of the necking type. For the other two materials, the post-critical solution associated with the only bifurcation point must necessarily be of the bulging type.

We also observe that the maximum or minimum of y occurs when y ′ = 0 and it occurs at Z = 0. It then follows from (4.14) that the maximum/minimum is given by

y max/min = -ω ′ (r cr )(r ∞ -r cr )/γ(r cr ). The relation r 0 ≈ r ∞ + y max/min = r ∞ -ω ′ (r cr )(r ∞ -r cr )/γ(r cr )
then provides a consistency check on our numerical values of r 0 . We conclude this section by remarking that the bifurcation problem addressed in this section corresponds to a bifurcation with zero mode number. Such bifurcations are usually associated with infinite domains but the corresponding spatially localized solutions are also valid for a finite domain since the finite boundaries (such as the ends of the tube) are not felt by the localized solutions. For such bifurcation problems, linear bifurcation analysis would incorrectly predict the eigenfunction as a constant; the correct shape of the eigenfunction can only be deduced from a weakly nonlinear, multiple-scale analysis. Such an analysis would usually yield an amplitude equation of the form

A XX + c 1 A + c 2 A 2 = 0, (4.24)

A c c e p t e d m a n u s c r i p t

where X is a far-distance variable, A is the amplitude function and c 1 and c 2 are constants; see Fu [29, p.362]. Multiplying (4.24) by A X and then integrating, we would obtain

(A X ) 2 + c 1 A 2 + 2 3 c 2 A 3 = const, (4.25)
which is seen to be of the same form as (4.13) derived earlier. It is due to the simplicity of the current problem that we were able to derive (4.13) without having to conduct a multiple-scale analysis in the standard manner.

Numerical solutions

Following Kyriakides and Chang [START_REF] Kyriakides | The initiation and propagation of a localized instability in an inflated elastic tube[END_REF], we may write down the following system of differential equations for λ 1 , λ 2 , φ:

λ ′ 1 = λ 2 sin φ, λ ′ 2 = Ŵ1 -λ 2 Ŵ12 Ŵ22 sin φ, (5.1)
φ ′ = Ŵ1 Ŵ2 cos φ - Ŵ ∞ 1 r ∞ Ŵ2 λ 1 λ 2 ,
where φ is defined in Fig. 1. These equations are solved subject to the initial conditions

λ 1 (0) = r(0), λ 2 (0) = z ′ (0), φ(0) = 0. ( 5.2)
We could also use λ 1 , λ 2 , r ′ as dependent variables (Yang and Feng [START_REF] Yang | On Axisymmetrical Deformations of Nonlinear Membranes[END_REF]), in which case we solve

λ ′ 1 = r ′ , λ ′ 2 = Ŵ1 -λ 2 Ŵ12 Ŵ22 r ′ λ 2 , (5.3)
(r ′ ) ′ = Ŵ1 λ 2 Ŵ2 + Ŵ1 -λ 2 Ŵ12 Ŵ22 - Ŵ1 λ 2 Ŵ2 r ′ λ 2 2 - Ŵ ∞ 1 r ∞ Ŵ2 λ 1 λ 2 2 1 - r ′ λ 2 2
subject to the conditions

λ 1 (0) = r(0), λ 2 (0) = z ′ (0), r ′ (0) = 0. ( 5.4)
For the Varga material, (2.7) and (2.5) reduce to

√ r ′2 + z ′2 = 2r ∞ r(r 2 ∞ + 2 -rr ∞ ) , (5.5)
and respectively. We then obtain from (5.5) a single first-order differential equation

z ′ √ r ′2 + z ′2 1 - 1 r(r ′2 + z ′2 ) - P r 2 2H = 1 - 1 r ∞ - P r 2 ∞ 2H , (5.6)
r ′ = ± 4r 2 ∞ r 2 (r 2 ∞ + 2 -rr ∞ ) 2 -z ′2 , (5.7)
where z ′ can be expressed explicitly in terms of r with the use of (5.5) and (5.6). We note that this simplification is possible because r ′ appears in (5.6) only through √ r ′2 + z ′2 which can be eliminated with the aid of (5.5).

We cannot solve (5.7) subject to the condition r(0) = r 0 directly since the latter is a fixed point of (5.7), but locking to the fixed point can be avoided by slightly perturbing the initial condition. For instance, we may use r(0) = r 0 + 0.0001 as the initial condition and the resulting solution provides a good enough approximation to the actual solitary wave type solution. This provides a useful check on our numerical solutions.

In Fig. 5(a,b), we have shown two post-critical solutions for the Gent material corresponding to r ∞ = 1.6 and r ∞ = 9.7, which are close to the two bifurcation values 1.6907 and 9.6523, respectively. We see good agreement between the numerical solutions and the asymptotic solution given by (4.21).

Earlier in Section 3 we have already indicated that not every initial condition at Z = 0 will correspond to a localized solution. The criterion (4.22), together with the plots in Figures 2(b), 3(b) and 4(b), enables us to deduce whether a solution will be localized when r ∞ is sufficiently close to a bifurcation value. For other values of r ∞ , we have to resort to numerical integrations.

We first discuss the Varga material and refer to Fig. 2. The criterion (4.22) indicates that the bifurcated solution is a localized bulge for r ∞ close to and less than r cr . As r ∞ is reduced from r cr , numerical integration of (5.1) shows that the solution remains localized only for r * < r ∞ < r cr , where, as pointed out by Haughton [START_REF] Haughton | Elastic Membranes, chapter 7 in Nonlinear Elasticity: Theory and Applications[END_REF], the lower bound r * is the value of r ∞ at which r(0), calculated according to (3.10), makes the right hand side of (5.5) blow up, that is r(0) = 2(1 + r 2 ∞ )/r 3 ∞ = r ∞ + 2/r ∞ from which we obtain r * = 2 1/4 ≈ 1.1892 (Haughton's value of r * ≈ 1.22 seems to be in error). For values of r ∞ outside the interval [r * , r cr ] the solution is either unbounded or periodic.

For the Gent material with J m = 97.2, our numerical calculation yields the results shown in Fig. 3, where again we have used solid lines to indicate association with localized solutions and dashed lines with either unbounded or periodic solutions. We find that the solid line segment originating from the first bifurcation point terminates at the nearest turning point. However, contrary to our expectation, the other solid line segment originating from the second bifurcation point terminates at r ∞ = 9.8046, before the nearest turning point at r ∞ = 9.91015 is reached. We have also carried out calculations for J m = 30 and have found qualitatively similar results.

Finally, for the Ogden material the amplitude diagram of which is given by Fig. 4(a), it is found that the end points that divide localized and unbounded/periodic solutions correspond to the turning point and the bifurcation value of r ∞ in Fig. 4(a). Thus, the bifurcated solution is localized only for 1.1525 < r ∞ < 1.6955. The upper and lower branches of the amplitude curve in Fig. 4(a) have asymptotic representations given by r 0 ∼ 1.23r 1.48 ∞ -0.39r 0.52 ∞ and r 0 ∼ 0, respectively. These asymptotic expressions are plotted in Fig. 4(a) in dashed lines in order to validate the numerical results.

Initiation pressure and further discussions

Stability analysis for the uniform inflation solution has been performed by Chen [START_REF] Chen | Stability and bifurcation of finite deformations of elastic cylindrical membranes -part I. stability analysis[END_REF]. For the case of pressure control and free ends, he found that for a stable solution the following conditions must be satisfied:

(λ 2 1 Ŵ11 -λ 1 Ŵ1 ) ≥ 0, Ŵ22 ≥ 0, Ŵ2 ≥ 0, λ 2 2 Ŵ22 (λ 2 1 Ŵ11 -λ 1 Ŵ1 ) -(λ 1 λ 2 Ŵ12 -λ 1 Ŵ1 ) 2 ≥ 0, (6.1)
whereas a strict inequality version of these conditions gives a sufficient condition for (local) stability.

We note that Ŵ22 > 0 can be interpreted as one of the conditions required by strong ellipticity for 2D deformations. The resultant axial force (in the membrane) required to maintain a constant λ 2 at infinity is F = 2πrhσ 2 = 2πH Ŵ2 , which has the same sign as Ŵ2 . It is known [START_REF] Corneliussen | Finite deformation of elastic membranes with application to the stability of an inflated and extended tube[END_REF][START_REF] Haughton | Bifurcation of inflated circular cylinders of elastic material under axial loading. I. Membrane theory for thin-walled tubes[END_REF] that, independently of the form of the strain energy function, F ≥ 0 if λ 1 λ 2 2 ≥ 1 and F ≤ 0 if λ 1 λ 2 2 ≤ 1. In particular, for λ 2 = 1 we must have F > 0 during inflation. Thus, we assume that the second and third inequalities in (6.1) are satisfied. The first inequality in (6.1) is then implied by the last inequality, and the latter, when evaluated at λ 1 = r ∞ , λ 2 = 1, is simply ω(r ∞ ) ≥ 0. Thus, the bifurcation condition ω(r ∞ ) = 0 marks We now argue that the first bifurcation pressure determined from ω(r ∞ ) = 0 can be used as the so-called initiation pressure that has featured in recent experimental studies. This is supported by the following observations: (i) it is at this pressure that the uniform inflation solution ceases to be stable; (ii) the near-critical bifurcated deformation is a bulge, of precisely the same form as observed experimentally; (iii) Figures 2(a), 3(a) and 4(a) show that the initial bulging is associated with a reduction in r ∞ from r cr , again just as observed experimentally. As further support, we have shown in Fig. 6 this pressure against the experimental data taken from figure 9 of Kyriakides and Chang [START_REF] Kyriakides | On the inflation of a long elastic tube in the presence of axial load[END_REF]. In Kyriakides and Chang's experiments, the tube was modeled as an Ogden material and the membrane tube was closed at one end at which an extra axial force F is also applied. For this case, our bifurcation condition ω(r ∞ ) = 0 would take the form

λ 2 2 Ŵ22 (λ 2 1 Ŵ11 -λ 1 Ŵ1 ) -(λ 1 λ 2 Ŵ12 -λ 1 Ŵ1 ) 2 = 0, (6.2)
and the scaled pressure and axial force would be given by

P R µH = 1 λ 1 λ 2 • Ŵ1 µ , F 2πµRH = Ŵ2 µ - λ 1 2λ 2 • Ŵ1 µ , (6.3)

A c c e p t e d m a n u s c r i p t

where the shear modulus µ is calculated according to 2µ = 3 n=1 µ n α n . For each specified λ 1 , we solve (6.2) to find the corresponding λ 2 , and hence the scaled pressure and axial force. We vary λ 1 in the interval [START_REF] Kydoniefs | Finite axisymmetric deformations of an initially cylindrical elastic membrane[END_REF][START_REF] Green | Large elastic deformations and non-linear continuum mechanics[END_REF] to generate the curve in Fig. 6. We see excellent agreement between our theoretical predications and the first four experimental results. It is not yet clear why the agreement becomes poorer for higher values of F , but we observe that the eight experimental results correspond to an axial stretch approximately equal to 1.12, 1.28, 1.49, 1.80, 2.23, 2.68, 3.13 and 3.49, respectively. Thus, the poor agreement corresponds to fairly high axial stretches.

Finally, we note that the methodology presented in the present paper can also be used to study nonlinear wave propagation in fluid-filled distensible tubes. Our amplitude equation (4.13) is in fact an integral of the well-known Kortewed-de Vries (KdV) equation u t -6uu x + u xxx = 0 (6.4) when u is assumed to take the form u = (γ/2)y(x -ct), c = ω ′ (r cr )ǫr 1 . The KdV equation is the prototypal evolution equation for waves where weak nonlinearity and weak dispersion operate simultaneously, and is known to have solitary-wave solutions. In a recent study by Epstein and Johnston [START_REF] Epstein | On the exact speed and amplitude of solitary waves in fluid-filled elastic tubes[END_REF], the authors found that no solitary waves could propagate if the fluid was initially stationary prior to wave propagation. This is in contrast with previous studies by Demiray [START_REF] Demiray | Solitary waves in prestressed elastic tubes[END_REF], Demiray and Dost [START_REF] Demiray | Axial and transverse solitary waves in prestressed thin elastic tubes[END_REF], and Erbay et al [START_REF] Erbay | Wave propagation in fluid-filled non-linear viscoelastic tubes[END_REF] where (6.4) was derived and the contrary was found in the context of a small-amplitude and long-wavelength approximation. We note that Epstein and Johnston's [START_REF] Epstein | On the exact speed and amplitude of solitary waves in fluid-filled elastic tubes[END_REF] equations (6.4) and (6.6) can be cast in our form (4.1) and (4.2). Our bifurcation condition ω(r ∞ ) = 0 would be replaced by the leading-order dispersion relation ω(c) = 0 which determines the linear wave speed c corresponding to infinite wavelength. Our near-critical post-bifurcation analysis would then correspond to a small-amplitude and long-wavelength analysis, thus enabling a direct comparison with previous studies. This will be carried out in a separate study.

Figure 1 :

 1 Figure 1: Axisymmetric deformation of a thin-walled rubber tube

Figure 2 :

 2 Figure 2: Amplitude r 0 -r ∞ and ω(r ∞ ) as functions of r ∞ for the Varga material. The dotted lines in (a) indicate the fact that the corresponding solution is either periodic or unbounded.

Figure 3 :

 3 Figure 3: Amplitude r 0 -r ∞ and ω(r ∞ ) as functions of r ∞ for the Gent material with J m = 97.2. The dotted lines in (a) indicate the fact that the corresponding solution is either periodic or unbounded.

Figure 4 :

 4 Figure 4: Amplitude r 0 -r ∞ and ω(r ∞ ) as functions of r ∞ for the Ogden material. The dotted lines in (a) indicate the fact that the corresponding solution is either periodic or unbounded.

Figure 5 :

 5 Figure 5: Near-critical bulging/necking solutions corresponding to the Gent material with J m = 97.2. (a) r ∞ = 1.6; (b) r ∞ = 9.7. Solid lines: numerical solution; dotted lines: asymptotic solution given by (4.21).

A c c e p t e d m a n u s c r i p t

Figure 6 :

 6 Figure6: Comparison of the bifurcation pressure (solid line) determined by the bifurcation condition (6.2) and the initiation pressure (circles) determined experimentally by Kyriakides and Chang[START_REF] Kyriakides | On the inflation of a long elastic tube in the presence of axial load[END_REF] and given in their figure9.

A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t A c c e p t e d m a n u s c r i p t

Acknowledgements

This work is supported by a joint grant from BBSRC and EPSRC (BB/D014786/1) under their Stem Cell Science and Engineering Initiative. We thank Professor Andrej Ilichev for some stimulating discussions.


Wolfgang Steiner
email: wsteiner@fh-wels.at

A Finite Element Based Stability Test for Equilibria of Flexible Structures in Circular Orbits

Keywords: Stability, Nonlinear Mechanical System, Relative Equilibrium, Finite Elements, Amended Potential, Space Elevator, Reduced Energy Momentum Method

Large flexible structures could be one of the most challenging space technologies of the future. For example, the well known space elevator may one day guide astronauts along a 36,000 km long tether from the earth's surface to a geostationary orbit. In this paper, we imagine an arbitrary space structure in a stationary circular motion around a celestial body.

Usually such systems require two problems to be solved: the strength of the materials and an overall stability analysis. The latter task is a non-trivial one since in orbiting systems not the angular rate, but the angular momentum must be kept constant when applying, for example, Dirichlet's criterion. In particular, this fact becomes important if the system's dimension is in the magnitude of the orbital radius. Therefore, this paper will focus on a general applicable method to determine the stability of flexible earth orbiting structures.

The analysis process is based on the reduced energy momentum method presented by J. C. Simo, T. A. Posbergh and J. E. Marsden, which has been customized for use in systems with cyclic coordinates by the author. In the presented approach the second variation of an amended potential is derived from a standard finite element analysis. An efficient stability test is introduced which limits the computational effort to the evaluation of a few eigenvalues and eigenvectors of the system's tangent stiffness matrix. Finally, two examples are discussed to demonstrate the application of the method.

A c c e p t e d m a n u s c r i p t 1 Introduction

Many technical systems operate in a so-called relative equilibrium configuration, which is characterized by the fact that the only dynamic process is a uniform motion in a symmetry group like a rotation about an axis at a constant angular rate. For example, consider a rotating turbine or a spacecraft in a stationary circular motion around a celestial body. To find a relative equilibirium for such systems one usually introduces a properly rotating reference frame in which the field of centrifugal forces acts as an external load. If the system consists of flexible bodies finite element programs can then be used to compute deformations and stresses by carrying out a standard static analysis.

However, for technical applications it is also necessary to analyze the stability of a relative equilibrium. For example, it is known that large space structures may become unstable, if the system's dimension is in the magnitude of the orbital radius. Here one must always take care about the fact that for a perturbation motion of a relative equilibrium not the velocity of the group parameter but the correspoding momentum is conserved. Therefore, in the case of a rotating system the angular velocity must not be treated as a constant quantity.

In this paper a general stability test for flexible orbiting structures is presented, taking into accont the correct conservation of the angular momentum. The test employs the method of the amended potential presented in [(5)] and [(4)] which has also been derived in [ [START_REF] Steiner | Stability Analysis of Relative Equilibria of Mechanical Systems with Cyclic Coordinates: A Direct Aproach[END_REF]] from an engineer's point of view. It can be applied as a post-processing step at the end of a standard finite element analysis and is also suitable for systems with many degrees of freedom since the computational effort is minimized by a proper decomposition of the second derivatives of the amended potential. Thus, only the eigenvalues and eigenvectors of the tangent stiffness matrix and the load vector of the centrifugal forces for the relative equilibrium are required as input for the test.

To demonstrate how the test works, the stability of a so-called space elevator is analyzed. One day this structure could enable astronauts to climb up a 36,000 km long tether connecting the earth's surface with a geostationary orbit. However, it turns out that an orbiting tether becomes unstable beyond a critical length [( 7), [START_REF] Steindl | Is the Sky-Hook Configuration Stable?[END_REF]].

Stability of relative equilibria

We briefly recall the results of the stability analysis process described in [ [START_REF] Steiner | Stability Analysis of Relative Equilibria of Mechanical Systems with Cyclic Coordinates: A Direct Aproach[END_REF]].

A relative equilibrium is closely related to the existence of cyclic coordinates,

A c c e p t e d m a n u s c r i p t

i.e. generalized coordinates that do not appear explicitly in the Lagrangian. In a relative equilibrium these coordinates vary uniformly like the rotation angle of a rotating system. Conversely, the non-cyclic coordinates remain stationary during the motion.

Let q c and q n be the vectors of the cyclic and the non-cyclic coordinates. The kinetic energy of the system is then given by

T = 1 2 qT M(q n ) q = 1 2 qT c , qT n    M cc M cn M nc M nn       qc qn    (1)
introducing M(q n ) as the mass matrix. For the stability analysis only the submatrix M cc , known as locked inertia tensor, is needed. It is derived from the mechanical model by setting the velocities of the non-cyclic coordinates qn to zero and qc = ω = const. In this case the kinetic energy reads

T c = 1 2 ω T M cc (q n )ω . (2)
Now, since qc is a cyclic coordinate, the corresponding generalized momentum

p c = ∂T ∂ qc = M cc (q n )ω (3)
is a first integral of the system and will therefore be conserved. This must be taken into account when disturbing a relative equilibrium in order to check its stability.

In Steiner [ [START_REF] Steiner | Stability Analysis of Relative Equilibria of Mechanical Systems with Cyclic Coordinates: A Direct Aproach[END_REF]] it is shown that relative equilibria and the correct stability test can be derived from the amendend potential

V * (q n ) = 1 2 p T c M -1 cc p c + V (q n ) (4)
where V (q n ) is the physical potential energy (depending on q n only). A set of nonlinear equations for the equilibrium configuration of the non-cyclic coordinates is obtained from

∂V * ∂q n = 0 at q n = q n,0 . (5)
Furthermore, the equilibrium is stable if the matrix of second derivatives of

V * A ij = ∂ 2 V * ∂q i n ∂q j n (6)
is positive definite at q n = q n,0 . We point out, that for all derivatives in ( 5) and ( 6) p c must be kept constant. Therefore, it is strictly forbidden to insert (3) into (4) before carrying out these operations.

A c c e p t e d m a n u s c r i p t 3 Modelling an orbiting flexible structure with finite elements

Let us consider now an arbitrary flexible structure rotating around a celestial body. Usually one introduces a rotating reference frame as depicted in fig. 1 to describe (relative) equilibrium configurations of the structure, since no timedepending coordinates will appear in this case. We can define the reference frame by a selected point that is constrained to the y-z-plane (x = 0). In fig. 1 this point is labelled "Node 1". Additionally, the z-Axis is fixed in the inertial frame.

Next we introduce a finite element discretization for the flexible structure.

Then the equilibrium configuration is described by a vector u ∈ R N of nodal degrees of freedom representing translations and rotations of the nodes with respect to an unstressed reference configuration. If we only include the unconstrained degrees of freedom in u, then this vector represents the non-cyclic coordinates of the system, i. e. q n = u. Note that the translation of node 1 along the x-axis is not an element of u since this degree of freedom is locked in order to define the reference frame. Now, to describe the configuration of the structure in an inertial frame we have to consider the rotation angle ϕ of the reference frame as an additional variable. In a relative equilibrium ϕ is a uniformly varying quantity, i. e.

φ = ω = const. (7)
Therefore, ϕ is the cyclic variable of the system.

Moreover, since q c = ϕ is the only cyclic variable, the locked inertia tensor M cc is a scalar M cc which depends on the non-cyclic coordinates u. As described in (2) it can be derived from the kinetic energy for u = 0. Since the system rotates around the z-axis we obtain

T c = 1 2 m v 2 dm = 1 2 m (ωr z ) 2 dm = 1 2 ω 2 m (x 2 + y 2 ) dm (8)
where x and y denote the x-and y-coordinate of a mass element dm in the deformed configuration. Comparing (8) with (2) we get

M cc = m (x 2 + y 2 ) dm (9)
where x and y can be expressed by u and the ansatz-functions introduced for the discretization process. Physically speaking, M cc = M cc (u) is the moment of inertia about the z-axis for a configuration u.

To describe the dynamic behaviour of our system the potential energy V (u) must be known. It consists of two parts, the strain energy U and the potential

A c c e p t e d m a n u s c r i p t

of the external load V G . The strain energy is the key functional used in finite element solvers for the derivation of "internal nodal forces". The external load of space vehicles is usually dominated by gravitational forces. Let κ be the gravitational parameter of the celestial body, then

V G = - m κ r dm = - m κ √ x 2 + y 2 + z 2 dm . (10)
Again, for numerical solutions x, y and z must be discretized by proper ansatzfunctions. Since not all software products allow the implementation of gravitational forces, it may become necessary to write user subroutines to take into account this type of loading, see section 6.

From the discretized functionals M cc and V = U + V G the amended potential (4) can be derived:

V * (u) = V (u) + 1 2 p 2 c M cc (u) . (11)
Here, p c is the conserved angular momentum

p c = M cc (u)ω (12)
which must be treated as a constant parameter in (11). Note, that all matrix products in (3) and ( 4) are reduced to simple one-dimensional operations in (11) and (12) since ϕ is the only cyclic variable.

Using the amended potential, the equations, which determine the relative equilibria u = u 0 , read

∂V * ∂u = ∂V ∂u - 1 2 p 2 c M 2 cc ∂M cc ∂u = 0 . (13)
In these equations the expression

f c = 1 2 p 2 c M 2 cc ∂M cc ∂u = 1 2 ω 2 ∂M cc ∂u (14)
is the load vector of centrifugal forces acting on the structure in the rotating reference frame. Most finite element solvers allow this type of load to be applied for a given ω. It is also easy to extract this vector from the FE-solver. For example, after computing the equilibrium in the first analysis step one can switch off the centrifugal loading while suppressing further nodal displacements in a subsequent step. The vector of nodal reaction forces (and moments) will then contain the elements of f c .

Proceeding from the amended potential we are now ready to analyze the stability of the relative equilibrium. To build up the matrix of second derivatives of V * the left side of (13) must be differentiated once more. Keeping p c constant

A c c e p t e d m a n u s c r i p t

we obtain

A ij = ∂ 2 V * ∂u i ∂u j = ∂ 2 V ∂u i ∂u j - 1 2 p 2 c M 2 cc ∂ 2 M cc ∂u i ∂u j + p 2 c M 3 cc ∂M cc ∂u i ∂M cc ∂u j . ( 15
)
Note that the last term would disappear if ω instead of p c were kept constant. However, we may substitute p c = M cc ω in (15) to describe A ij in terms of ω:

A ij = ∂ 2 V ∂u i ∂u j - 1 2 ω 2 ∂ 2 M cc ∂u i ∂u j + ω 2 M cc ∂M cc ∂u i ∂M cc ∂u j . ( 16
)
If the system is stable the matrix A = {A ij } must be positive definite for u = u 0 .

Stability analysis

It is the main goal of this paper to derive a simple stability test for flexible orbiting structures from ( 16), which uses only data available from a standard finite element analysis. Therefore, we note that A ij seperates into two parts.

First, the matrix

K ij = ∂ 2 V ∂u i ∂u j - 1 2 ω 2 ∂ 2 M cc ∂u i ∂u j ( 17
)
is the tangent stiffness matrix of the system for a constant rotation of the reference frame. It includes all structural stiffness effects as well as load stiffness contributions due to centrifugal and gravitational forces. The FE-software also uses this matrix to solve for equilibrium configurations by means of the Newton-Raphson method. It can be exported from most FE-programs. However, we shall see that this is not necessary.

The latter part of ( 16)

ω 2 M cc ∂M cc ∂u i ∂M cc ∂u j = ω √ M cc ∂M cc ∂u i ω √ M cc ∂M cc ∂u j (18)
is the tensor product of the vector

a = ω √ M cc ∂M cc ∂u = ω √ M cc 2 ω 2 f c = 4 M cc ω 2 f c = 2 T c f c ( 19
)
where we have used ( 14) and (2), i. e.

T c = M cc ω 2 /2.
With (17), ( 18) and ( 19) we may now rewrite (16) in the following way:

A ij = K ij + a i a j ⇐⇒ A = K + a a T (20)

A c c e p t e d m a n u s c r i p t

Here, the matrix notation aa T = a⊗a was used to express the tensor product.

A is positive definite, if the quadratic form

x T Ax = x T K + a a T x (21)
is positive for an arbitrary vector x ∈ R N and zero only for x = 0. In this case, the system is stable. Now, to show that A is positive definite we could simply compute the lowest eigenvalue of A and see if it is positive. However, A contains a tensor product and is therefore a dense matrix. For large FE-models the eigenvalue extraction of A may cause difficulty. On the other hand, the tangent stiffness matrix K is a sparse matrix for which eigenvectors and eigenvalues can be extracted efficiently.

Let λ 1 , . . . λ N be the eigenvalues of K sorted by its magnitude, i. e. λ 1 < λ 2 < • • • < λ N , and let v 1 , . . . v N be the corresponding eigenvectors of K.

(For example, these objects can be obtained from FE-solvers by performing a natural frequency extraction, if the mass matrix is set to the unity matrix. Note that the base state of this extraction is the relative equilibrium.) Since K is symmetric, see ( 17), the eigenvalues are real numbers and the eigenvectors form an orthonormal basis of R N . Thus,

v T i v i = 1 and v T i v j = 0 for i = j . (22)
In other words, the N × N matrix

V = (v 1 , v 2 , . . . v N ) (23)
is orthogonal, i. e. V T = V -1 . Furthermore, the stiffness matrix diagonalizes under the following operation:

V T KV = D = diag(λ i ) . (24)
The transformation x = Vy simply describes a change of basis for vectors in R N . We can use it to rewrite the quadratic form (21):

F (y) = y T V T K + a a T Vy = y T V T KV + V T a a T V y
Inserting (24) and introducing the vector b = V T a, which has the components

b i = v T i a , ( 25
)
this expression is reduced to

F (y) = y T D + b b T y . (26)

A c c e p t e d m a n u s c r i p t

D is a diagonal matrix containing the eigenvalues of K and the components of b are the projections of the vector a onto the eigenvectors of K. So our quadratic form is defined by the matrix:

B = D + b b T =           λ 1 + b 2 1 b 1 b 2 • • • b 1 b N b 2 b 1 λ 2 + b 2 2 • • • b 2 b N . . . . . . . . . . . . b N b 1 b N b 2 • • • λ N + b 2 N           (27)
Now, B (and therefore A) is positive definite, if and only if the determinants of all k × k submatrices,

B k =               λ 1 + b 2 1 b 1 b 2 b 1 b 3 • • • b 1 b k b 2 b 1 λ 2 + b 2 2 b 2 b 3 • • • b 2 b k b 3 b 1 b 3 b 2 λ 3 + b 2 3 • • • b 3 b k . . . . . . . . . . . . . . . b k b 1 b k b 2 b k b 3 • • • λ k + b 2 k              
are positive (leading principal minor criterion). In the appendix it is shown, that

det B k = λ 1 λ 2 . . . λ k 1 + b 2 1 λ 1 + b 2 2 λ 2 + • • • + b 2 k λ k . ( 28
)
Hence, the relative equilibrium computed from (13) is stable, if det B k > 0 for k = 1, . . . N.

Fast stability tests

A general stability test for a rotating flexible structure might still be time consuming since all eigenvectors and eigenvalues of the tangent stiffness matrix are needed to show that det B N > 0. In particular, this could cause difficulties if the structure has many degrees of freedom. In such a case faster tests which enable the stability to be checked without computing det B k from (28) for all k = 1, . . . N need to be found.

Basically, the sign of det B k is determined by the signs of the eigenvalues

λ 1 < λ 2 < • • • < λ N .
Thus, we consider the following cases:

A c c e p t e d m a n u s c r i p t

Case 1: λ 1 > 0

Stability is granted, if λ 1 is positive because then all det B k are positive. For this check, only the lowest eigenvalue of K must be computed.

Case 2: λ 1 < 0, λ 2 < 0

In this case at least two eigenvalues of K are negative. Suppose the system is stable, then

det B 1 = λ 1 1 + b 2 1 λ 1 > 0
and since λ 1 < 0 the bracket must also be negative. Now look at

det B 2 = λ 1 λ 2 1 + b 2 1 λ 1 + b 2 2 λ 2 .
Since λ 2 < 0 the product λ 1 λ 2 is positive and the sign is determined by the bracket. But as seen above 1 + b 2 1 /λ 1 < 0. Because also b 2 2 /λ 2 < 0 we obtain det B 2 < 0. Thus, the system is unstable. For this check, only the first two eigenvalues of K must be computed.

Case 3: λ 1 < 0, λ 2 > 0
The most complex situation arises, if only one eigenvalue is negative. In this case the factor λ 1 λ 2 . . . λ k in (28) is negative for all determinants. If det B k > 0, also the brackets must be negative:

1 + b 2 1 λ 1 + b 2 2 λ 2 + • • • + b 2 k λ k < 0 for k = 1, . . . N.
On the left side only the term b 2 1 /λ 1 is negative by assumption. So the worst case is obtained for k = N, where all terms b 2 i /λ i are added. We may conclude that the system is stable only if

1 + b 2 1 λ 1 + b 2 2 λ 2 + • • • + b 2 N λ N = 1 + N i=1 b 2 i λ i < 0 . (29)
However, one would prefer an iterative stability test to avoid the computation of all eigenvalues and eigenvectors. In fact, if for any k < N

1 + k i=1 b 2 i λ i > 0 (30)

A c c e p t e d m a n u s c r i p t

the system is unstable and the stability test can be stopped at that point. On the other hand, if the sum in ( 29) is truncated at i = k we can also find an upper limit for the missing terms:

N i=k+1 b 2 i λ i ≤ 1 λ k+1 |b| 2 - k i=1 b 2 i since λ k+1 < λ k+2 < • • • < λ N and b 2 k+1 ≤ |b| 2 -k i=1 b 2 i . Thus, the system is stable if for any k < N 1 + k i=1 b 2 i λ i + 1 λ k+1 |b| 2 - k i=1 b 2 i < 0 (31)
For the tests (30) and (31) only the first k + 1 eigenvalues and the first k eigenvectors must be computed. Note that

|b| 2 = |a| 2
since b is just the decompostion of a in the orthonormal basis v i .

The complete stability test including all cases is depicted in the structogram of fig. 2.

Remarks to the finite element method

The complete stability analysis can be carried out with standard FE-software packages. We describe an ABAQUS based approach in this section. In ABAQUS gravitational and centrifugal forces should be included via additional userelements in order to introduce their correct contribution to the residuum and to the stiffness matrix. Such a user-element is added for every structural finite element.

Basically, the Potential V G from (10) and the locked inertia tensor M cc from (9) must be expressed by the vector u of nodal degrees of freedom. This is usually done by subdividing the structure into finite elements. Let u (e) denote the vector of nodal degrees of freedom for element (e). Then the position r of a material point within the element, located at r 0 in the reference configuration, is given by r = Ψ (e) u (e) + r 0 (32)

where Ψ (e) is a linear operator mapping u (e) ∈ R N (e) onto r -r 0 ∈ R 3 . It is defined by the ansatz-functions used for the discretization process.

If a material point dm is described by a proper density ̺ and by an infinitesimal Ψ (e)T g( r) dΩ, where g( r) = ∂γ ∂ r = -̺κ r | r| 3 (34)

and the element load stiffness matrix

K (e) G,ij = ∂ 2 V (e) G ∂u (e) i ∂u (e) j is K (e) G = Ω (e)
Ψ (e)T G( r)Ψ (e) dΩ, where

G( r) = ∂ g ∂ r = - ̺κ | r| 3 3 r ⊗ r | r| 2 -I (35
) After numerical integration, (34) and ( 35) can be coded in the FORTRAN user-subroutine "UEL" provided by ABAQUS. In a similar way the locked inertia tensor may be considered. Altogether, the user elements yield the correct equilibrium equations and the tangent stiffness matrix K which is needed for the stability test.

Examples

The Space Elevator

To demonstrate an application of the stability test described above we first consider the "space elevator" depicted in fig. 3. For simplicity, the elevator is modeled as homogenous string without bending and torsional stiffness. The modelling space is three-dimensional. Thus, motions out of the orbital plane are taken into account.

For an inextensible string the radial equilibrium configuration, see fig. 3, can be derived directly from Newton's law. Suppose that the system rotates at a given angular rate ω around a celestial body exhibiting the gravitational parameter κ. Furthermore, let L and µ be length and mass per unit length of the string. Introducing R as the orbital radius of the string's mass center we obtain from Newton's second law

µLRω 2 = R+L/2 R-L/2 µκ r 2 dr = -µκ 1 R + L/2 - 1 R -L/2

A c c e p t e d m a n u s c r i p t

yielding

ω 2 = κ R 3 1 1 -1 4 L R 2 (36)
For the subsequent analysis L is increased from 0.1R to 1.8R. Simultaneously, ω is computed from (36), although the string is not completely inextensible in the finite element model. Furthermore, the parameters R, µ and κ are set to unity in order to normalize the problem. To find an appropriate value for the axial stiffness of the string the dimensionless parameter EAR/(µκ) can be computed for a realistic situation. For example, if the center of mass is in a geostationary orbit around the earth, then R = 42, 000 km and κ = 3.986 10 14 m 3 /s 2 . For a steel string Young's modulus is E = 2.1 10 11 N/m 2 and the density is given by ̺ = 7, 850 kg/m 3 . Thus, since µ = ̺A, the stiffness parameter would be

EA R µκ = ER ̺κ = 3.22 .
However, in our example the stiffness is increased to 1000.0 in order to approximate the inextensible case.

The result of the stability analysis is explained in fig. 4 and fig. 5. Fig. 4 shows the first two eigenvalues of the tangent stiffness matrix as a function of the string length. Obviously λ 1 is negative and λ 2 is positive in the considered interval. Therefore, case 3 of the stability test discussed in section 5 is relevant here.

As pointed out in section 5 the proof that the system is stable or unstable may require the computation of variable numbers k of eigenvalues. Fig. 5 shows that this number increases along with the string length. At L ≈ 1.0 thirteen eigenvalues of the tangent stiffness matrix must be evaluated. Beyond L ≈ 1.0 the system becomes unstable which is proven by the first eigenvalue λ 1 .

This result is in agreement with numerical simulations of the perturbation motions and also with a simplified analytical model of the space elevator. First, fig. 6 shows the perturbation motion of the lower endpoint of the string for a stable and for an unstable configuration. For that purpose a fully dynamic finite element simulation of the string was carried out. It can be observed, that for L = 0.946R (solid line) the point rotates around the earth whereas for L = 1.018R (dashed line) the system descends after a few revolutions. Since both configurations were chosen close to the stability boundary, we may conclude that the stability analysis yielded a reasonable result.

To verify the obtained result we may also derive an analytical stability condition for an inextensible string. According to [( 6)] it should be sufficient to consider radial perturbation motions only. Thus, we only have to vary the distance R of the string's mass center in the amended potential.

A c c e p t e d m a n u s c r i p t

In the radial configuration shown in fig. 3 the locked inertia tensor according to (9) is simply given by

M cc = m R 2 + L 2 12 ,
where m denotes the total mass of the string. The potential energy reads

V = - R+L/2 R-L/2 µκ r dr = -µκ ln R + L/2 R -L/2
resulting in the amended potential

V * (R) = p 2 2M cc (R) + V (R) ( 37
)
where p is the angular momentum of the string. By setting the first derivative of V * (R) to zero we obtain a relation between the angular momentum and the value of R = R 0 at the relative equilibrium:

p 2 = mκµ 36 L R 0 (L 2 + 12R 2 0 ) 2 4R 2 0 -L 2 . (38)
The stability condition for the inextensible string is then obtained from

d 2 V * dR 2 R=R 0 > 0 ( 39
)
which may be expressed in terms of R 0 by eliminating p using (38). Finally this procedure results in

4κµ L R 0 L 4 -48R 2 0 L 2 + 48R 4 0 (L 2 -4R 2 0 ) 2 (L 2 + 12R 2 0 ) > 0 .
This expression changes its sign if

L 4 -48R 2 0 L 2 + 48R 4 0 = 0 which happens at L = 2 6 - √ 33R 0 ≈ 1.01082R 0 .
This is in good agreement with the numerical stability test for the orbiting string shown in fig. 5.

The Orbiting Ring

The second example to be discussed here, is the flexible ring shown in fig. 7. For this structure also bending and torsional stiffness are taken into account.

A c c e p t e d m a n u s c r i p t

Again, the problem is treated in three-dimensional space. The basic parameters are selected as follows: Let E be Young's modulus, A the area and J the moment of inertia for bending of the ring's cross section. Assuming a thin walled circular cross section with radius a and thickness h, then A = 2πah and J = πa 3 h .

Furthermore, µ = ̺A denotes the mass per unit length and R the orbital radius of the reference configuration. It is convenient to formulate the problem by means of dimensionless quantities Ē, Ā, J, μ, κ and R. As for the orbiting string we set μ = 1, κ = 1 and R = 1. To find appropriate values for Ē, Ā and J we may consider dimensionless quantites such as

A R 2 = Ā R2 = Ā , J R 4 = J R4 = J ER 3 µκ = Ē R3 μκ = Ē
Therefore, the dimensionless axial and bending stiffness are given by

Ē Ā = EAR µκ = E ̺ R κ , Ē J = EJ µκR = E ̺ J ARκ .
Both expressions do not depend on the wall thickness h of the cross section. However, materials with very high specific Young's modulus E/̺ and a relatively large radius a of the cross section are necessary to avoid complete folding of the ring under the gravitational and centrifugal forces. We choose R = 42, 000 km and κ = 3.986 10 14 m 3 /s 2 representing a geostationary orbit around the earth, and set E/̺ = 1.9 × 10 10 Nm/kg, which is about 700 times higher than the respective value for steel, yielding Ē Ā = 2000. For the cross section radius we assume a = 0.01R = 420 km resulting in Ē J = 0.1. Although these are quite futuristic assumptions compared with today's techonology it is interesting to observe the stability behaviour of such structures.

To carry out the stability analysis the angular rate must be defined for the equilibrium. Since we are looking for a geostationary configuration we choose

ω 2 R 3 κ = ω2 R3 κ = ω2 = 1 .
Despite the vast dimensions of the structure the ring deforms significantly under the gravitational and the centrifugal load. Fig. 8 shows the ovalization for varying ring radius r = r/R. Note that the long axis is pointing towards the center of gravity. The stability analysis shows that the equilibrium is stable up to r ≈ 0.4. As for the orbiting string, λ 1 is negative and λ 2 is positive, see fig. 9. Thus, case 3 of the stability test has to be considered again. For small values of r only three eigenvalues are required to show that the system is stable.

A c c e p t e d m a n u s c r i p t

However, more eigenvalues must be computed in the vicinity of the stability bound (e. g. 80 at r ≈ 0.41). Beyond the stability bound the proof yields the result from the first eigenvalue. Fig. 10 shows the system at the stability boundary. For the undeformed reference configuration also the dimension of the cross section can be observed. Like in the string-example one may expect that loss of stability affects the overall orbital dynamics and that the system will descend.

Both examples described in this section demonstrate that in many cases the first few eigenvalues and eigenvectors of the tangent stiffness matrix K are sufficient to obtain information about the stability of an orbiting structure. Furthermore, the test can be integrated efficiently in a finite element analysis with common software packages.

APPENDIX: A formula for the determinant of B k

By mathematical induction we show that

det B k = det               λ 1 + b 2 1 b 1 b 2 b 1 b 3 • • • b 1 b k b 2 b 1 λ 2 + b 2 2 b 2 b 3 • • • b 2 b k b 3 b 1 b 3 b 2 λ 3 + b 2 3 • • • b 3 b k . . . . . . . . . . . . . . . b k b 1 b k b 2 b k b 3 • • • λ k + b 2 k               = λ 1 λ 2 . . . λ k 1 + b 2 1 λ 1 + b 2 2 λ 2 + • • • + b 2 k λ k (40)
where B k is the k × k-matrix introduced in section 4 for the formulation of the stability test.

B k = 1 b 2 k-1 det                  λ 1 + b 2 1 b 1 b 2 • • • b 1 b k-2 b 1 b k-1 b 1 b k b k-1 b 2 b 1 λ 2 + b 2 2 • • • b 2 b k-2 b 2 b k-1 b 2 b k b k-1 . . . . . . . . . . . . . . . . . . b k-2 b 1 b k-2 b 2 • • • λ k-2 + b 2 k-2 b k-2 b k-1 b k-2 b k b k-1
A c c e p t e d m a n u s c r i p t

det B k = 1 b 2 k-1 λ k b 2 k-1 + λ k-1 b 2 k det B k-1 + λ k-1 b k (-λ k-1 b k ) det B k-2 = b 2 k b 2 k-1 λ k-1 det B k-1 -λ k-1 det B k-2 + λ k det B k-1
If we now substitute the induction hypothesis for det B k-1 and det B k-2

det B k-1 = λ 1 λ 2 . . . λ k-1 1 + b 2 1 λ 1 + b 2 2 λ 2 + • • • + b 2 k-1 λ k-1 det B k-2 = λ 1 λ 2 . . . λ k-2 1 + b 2 1 λ 1 + b 2 2 λ 2 + • • • + b 2 k-2 λ k-2 we conclude that det B k = b 2 k b 2 k-1 λ k-1 λ 1 λ 2 . . . λ k-1 b 2 k-1 λ k-1 +λ 1 λ 2 . . . λ k 1 + b 2 1 λ 1 + b 2 2 λ 2 + • • • + b 2 k-1 λ k-1 = λ 1 λ 2 . . . λ k 1 + b 2 1 λ 1 + b 2 2 λ 2 + • • • + b 2 k λ k .
This is exactly the hypothesis for det B k which is therefore proved. Since we have supposed b k-1 = 0 for our proof we finally must consider the situation where some b i are zero. However, in this case the matrix B k can be rearranged by exchanging rows and columns so that all elements with non-zero b i are collected in the upper left corner of the matrix. Now let b1 , . . . bm denote the sequence of all non-zero b i appearing in B k and let λ1 , . . . λm be the corresponding sequence of λ i (m ≤ k). Then by exchanging rows and columns we obtain

det B k = det                      λ1 + b2 1 b1 b2 • • • b1 bk 0 • • • 0 b2 b1 λ2 + b2 2 • • • b2 bk 0 • • • 0 . . . . . . . . . . . . . . . . . . . . . bm b1 bm b2 • • • λm + b2 m 0 • • • 0 0 0 • • • 0 λm+1 • • • 0 . . . . . . . . . . . . . . . . . .

A c c e p t e d m a n u s c r i p t

Clearly, this determinant is given by det

B k = λm+1 . . . λk det           λ1 + b2 1 b1 b2 • • • b1 bk b2 b1 λ2 + b2 2 • • • b2 bk . . . . . . . . . . . . bm b1 bm b2 • • • λm + b2 m          
For the remaining determinant we may already use our hypothesis since per definition bi = 0. Introducing the sequence of zero b i 's by bm+1 , . . . bk = 0 we obtain det B k = λm+1 . . . λk λ1 . . . λm 1 + b2

1 λ1 + • • • + b2 m λm + b2 m+1 λm+1 + • • • + b2 k λk
In this expression the sequence is irrelevant and we may omit the bars yielding again the formula for det B k . Thus, our proof is complete. Fig. 6. Two orbits of the lower end point of the string from fig. 3 after slight perturbations of the radial equilibria. For L = 0.946R (solid line) the equilibrium is stable whereas for L = 1.018R (dashed line) the the equilibrium has become unstable. 28 Figure10

A c c e p t e d m a n u s c r i p t volumetric interval dΩ

  the gravity potential for a finite element reads V

First, ( 40 A c c e p t e d m a n u s c r i p t (

 40t ) is readily verified for k = 1 and k = 2. Then we have to show that 40) is true for any k if it is true for k -1 and k -2. This can be done by simplifying the determinant in the following way.Provided that b k-1 = 0, multiplying the last row and the last column of B k with b k-1 yields det

Fig. 1 .

 1 Fig. 1. A flexible orbiting structure in a rotating reference frame.

  19

Figure1 A c c e p t e d m a n u s c r i p t λ 1 Fig. 2 .

 12 Figure1

  20

Figure2Fig. 3 .

 3 Figure2

  21

Figure3

  Figure3

Fig. 4 .

 4 Fig. 4. First eigenvalues of the orbiting string for varying length L.

  22

Figure4

  Figure4

Fig. 5 .

 5 Fig. 5. Number of eigenvalues and eigenvectors of the orbiting string needed for the stability test.

  23

Figure5

  Figure5

Figure6Fig. 7 .

 7 Figure6

  25

Figure7Fig. 8 .

 8 Figure7

  26

Figure8 1 Fig. 9 .

 19 Figure8

  27

Figure9Fig. 10 .

 10 Figure9

b 2 b k-1 b 2 b k b k-1 . . . . . . . . .

Acknowledgements

I would like to thank my collegues Klaus Schiefermayr for the fruitful discussions of mathematical questions and Douglas Vaught for improving the linguistic quality of the paper. Furthermore, I thank the reviewers for the helpful comments on the discussion of the examples which were treated more extensivly after revision.