Zhu, Yunfei (2010) Nonlinear deformations of a thick-walled

hyperelastic tube under external pressure. PhD thesis.

http://theses.gla.ac.uk/1627/

Copyright and moral rights for this thesis are retained by the author

A copy can be downloaded for personal non-commercial research or

study, without prior permission or charge

This thesis cannot be reproduced or quoted extensively from without first

obtaining permission in writing from the Author

The content must not be changed in any way or sold commercially in any
format or medium without the formal permission of the Author

When referring to this work, full bibliographic details including the

author, title, awarding institution and date of the thesis must be given.

Glasgow Theses Service

http://theses.gla.ac.uk/
theses@gla.ac.uk
 Nonlinear deformations of a
thick-walled hyperelastic tube
under external pressure

by
Yunfei Zhu

A thesis submitted to the

Faculty of Information and Mathematical Sciences
at the University of Glasgow
for the degree of
Doctor of Philosophy

December 2009

c Yunfei Zhu 2009

°
Abstract

This research deals with several novel aspects of the nonlinear behaviour of thick-walled
cylindrical hyperelastic tubes under external pressure.
Initially, we consider bifurcation from a circular cylindrical deformed configuration of a
thick-walled circular cylindrical tube of incompressible isotropic elastic material subject to
combined axial loading and external pressure. In particular, we examine both axisymmet-
ric and asymmetric modes of bifurcation. The analysis is based on the three-dimensional
incremental equilibrium equations, which are derived and then solved numerically for a
specific material model using the Adams-Moulton method. We assess the effects of wall-
thickness and the ratio of length to (external) radius on the bifurcation behaviour.
The problem of the finite axisymmetric deformation of a thick-walled circular cylindri-
cal elastic tube subject to pressure on its external lateral boundaries and zero displacement
on its ends is formulated for an incompressible isotropic neo-Hookean material. The for-
mulation is fully nonlinear and can accommodate large strains and large displacements.
The governing system of nonlinear partial differential equations is derived and then solved
numerically using the C++ based object-oriented finite element library Libmesh. The
weighted residual-Galerkin method and the Newton-Krylov nonlinear solver are adopted
for solving the governing equations. Since the nonlinear problem is highly sensitive to
small changes in the numerical scheme, convergence was obtained only when the analyt-
ical Jacobian matrix was used. A Lagrangian mesh is used to discretize the governing
partial differential equations. Results are presented for different parameters, such as wall
thickness and aspect ratio, and comparison is made with the corresponding linear elas-
ticity formulation of the problem, the results of which agree with those of the nonlinear
formulation only for small external pressure. Not surprisingly, the nonlinear results depart
significantly from the linear ones for larger values of the pressure and when the strains in
the tube wall become large. Typical nonlinear characteristics exhibited are the “corner

i
 ii

bulging” of short tubes, and multiple modes of deformation for longer tubes.
Finally the general fully nonlinear governing equations in Lagrangian form for the three
dimensional large deformations of an elastic tube under external pressure are developed.
Acknowledgements

I would like to thank my supervisors, Professor Xiaoyu Luo and Professor Ray W. Ogden,
for their continued support and advice through out the period of this study. I feel extremely
fortunate to have been given the chance to work under them as a postgraduate and to
continue working with them as a research assistant. I am also grateful to the Overseas
Research Students Awards Scheme, and the Faculty of Information and Mathematical
Science, and the Department of Mathematics at Glasgow University for the generous
scholarships.
I also wish to thank Dr. David Haughton, University of Glasgow, and Professor Chris
Berttram, University of New South Wales, for their valuable advice on Chapter 4. I also
want to express my appreciation to Libmesh developers at the University of Texas and Dr.
Boyce Griffith at New York University for their help with using Libmesh.
My parents and elder sister deserve the deepest thanks for their never-ending support
during my education.
Finally, I would like to thank the cool kids in my office, such as Robert, Stephen,
Ehsan, Susan, Monica and Marjory for their help in the last three years.

iii
Contents

Abstract i

Acknowledgements iii

1 Introduction 1

2 Basic equations 11
2.1 Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 11
2.2 Stress theory and equilibrium . . . . . . . . . . . . . . . . . . . . . . . . . . 13
2.3 Constitutive law for a Cauchy elastic material . . . . . . . . . . . . . . . . . 14
2.4 Green elastic material . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.1 Stress-deformation relations in terms of invariants . . . . . . . . . . 17
2.5 Internal constraints . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.6 Example strain-energy functions for isotropic elastic material . . . . . . . . 19
2.7 Incremental deformations . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

3 Finite element method, libmesh library 22

3.1 Finite element nonlinear analysis in solid mechanics . . . . . . . . . . . . . 22
3.1.1 Introduction . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.1.2 Procedure for finite element solution with libmesh . . . . . . . . . . 23
3.1.3 Stress recovery . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.2 Libmesh library . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.2.1 Object-Oriented Scientific Computing . . . . . . . . . . . . . . . . . 28
3.2.2 libmesh: Adaptive mesh refinement scheme . . . . . . . . . . . . . . 29
3.2.3 Data structure of libmesh . . . . . . . . . . . . . . . . . . . . . . . . 30

iv
CONTENTS v

4 Bifurcation 34
4.1 The elastic constitutive law and strain-energy function . . . . . . . . . . . . 34
4.2 The circular cylindrical configuration . . . . . . . . . . . . . . . . . . . . . . 35
4.3 Incremental equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 37
4.4 Asymmetric bifurcations and numerical methods . . . . . . . . . . . . . . . 40
4.5 Numerical results and discussion . . . . . . . . . . . . . . . . . . . . . . . . 44
4.5.1 Equilibrium pressure curves . . . . . . . . . . . . . . . . . . . . . . . 44
4.5.2 Axisymmetric bifurcation . . . . . . . . . . . . . . . . . . . . . . . . 45
4.5.3 Asymmetric bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.5.4 Discussion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 56
4.6 Conclusion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58

5 Nonlinear axisymmetric deformations 60

5.1 Basic equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.1.1 Deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
5.1.2 Material properties and equilibrium . . . . . . . . . . . . . . . . . . 61
5.2 Linear and nonlinear equations . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.2.1 Radially symmetric case . . . . . . . . . . . . . . . . . . . . . . . . . 63
5.2.2 The linear case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 65
5.2.3 The nonlinear case . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66
5.3 Finite element algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.3.1 Discretization . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.3.2 Newton’s method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68
5.3.3 Detailed discretizing integrations . . . . . . . . . . . . . . . . . . . . 68
5.4 Numerical algorithms . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 70
5.5 Numerical results . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
5.5.1 Thick-walled short tubes: A/B = 0.5 and L/B = 1 . . . . . . . . . . 71
5.5.2 Thick-walled longer tubes: A/B = 0.5 and L/B = 5 . . . . . . . . . 87
5.5.3 Thinner and longer tubes: A/B = 0.8 and L/B = 5 . . . . . . . . . 90
5.6 Discussion and conclusions . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

6 Three-dimensional large deformations 94

6.1 Nonlinear case: cylindrical polar coordinates . . . . . . . . . . . . . . . . . . 94
6.1.1 Deformation gradient . . . . . . . . . . . . . . . . . . . . . . . . . . 94
CONTENTS vi

6.1.2 Nominal stress and Cauchy stress . . . . . . . . . . . . . . . . . . . . 96

6.1.3 Equilibrium equations . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.2 Nonlinear case: cartesian coordinates . . . . . . . . . . . . . . . . . . . . . . 97
6.2.1 Deformation gradient . . . . . . . . . . . . . . . . . . . . . . . . . . 97
6.2.2 Nominal stress and Cauchy stress . . . . . . . . . . . . . . . . . . . . 98
6.2.3 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 98
6.2.4 Unit outward normal to the outer surface . . . . . . . . . . . . . . . 99
6.3 Linear case . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99
6.3.1 Deformation gradient and Cauchy stress . . . . . . . . . . . . . . . . 100
6.3.2 Governing equations . . . . . . . . . . . . . . . . . . . . . . . . . . . 100

7 Conclusions 102

A Appendix A 105
A.1 Integration of the equilibrium equations by parts . . . . . . . . . . . . . . . 105

B Appendix B 109
B.1 Mesh files . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 109

References 113
Chapter 1

Introduction

In this thesis, we focus on modelling and simulating the collapse of a cylindrical tube
under external pressure. The collapse of a circular tube is of interest not only in many
engineering applications, such as in the design of undersea or underground pipelines and
pressure vessels, particularly submarine structures, but it is also of great interest in the
biomechanics context.
The main motivation for this research comes from the investigation of the behaviour
of physiological conduits within our body subject to external pressure, such as veins and
arteries conveying blood flow and large bronchi conducting air into the lungs. Due to
the high flexibility of the soft biological tissue, these conduits may collapse under certain
conditions of external and internal fluid pressure. For example, intramyocardial arteries
collapse during systole [29]. More generally, the collapse of cardiovascular vessels plays an
important role in the delivery of blood to other organs [30], [63]. A long cylindrical elas-
tic tube when subjected to a transmural (internal minus external) pressure may collapse
into a two-lobed configuration, as illustrated in Fig.1.1. In elastic buckling the collapse is
usually sudden and catastrophic, which is a process involving some nonlinear dynamical
deformations of great complexity. Depending on geometry, material properties, the pres-
sure and boundary conditions, the tubes may collapse differently. In essence, these kinds
of problems involve two physical systems interacting with each other, which are the elastic
wall of the conduit and the biological fluid inside or around the conduit. Such systems are
also known as coupled and such coupling may be weak or strong depending on the degree
of interaction.
A general definition of coupled systems may be given as [80]

1
CHAPTER 1. INTRODUCTION 2

Figure 1.1: Collapsed tube (from paper by Marzo et al. [52]).

Definition 1.0.1 Coupled systems and formulations are those applicable to multiple do-
mains and dependent variables which usually but not always describe the different physical
phenomena and in which
(a) neither domain can be solved while separated from the other;
(b) neither set of dependent variables can be explicitly eliminated at the differential
equation level.

The topic of flow through collapsible tubes, an obvious coupled problem, has been
studied for several decades. This subject has been reviewed briefly by Kamm and Pedley
[42]. Experiments [14], [12] on a Starling resistor prototype of the system have presented
a rich dynamic behaviour, with various types of self-excited oscillations. One- or two-
dimensional theoretical models has been established by Pedley [61], Luo [49], [50], [51],
Cai [17] and Jensen [41] and much computational work has been carried out to reveal the
mechanisms of such oscillations. Luo and Pedley [49] studied steady flow in a 2-D channel
with one plane rigid wall and the other wall replaced by an elastic segment, which is
treated as a elastic membrane, as illustrated in Fig.1.2. Following their previous work [49]
on steady flow in a two-dimensional fluid-membrane model of the collapsible tube, Luo
and Pedley [50] investigated the instability of the steady solution by developing a time-
dependent simulation of the coupled flow-membrane problem. These studies provided
a detailed picture of the fluid and solid mechanics involved in the large-amplitude self-
excited oscillations in this simplified system and have shown rich dynamical behaviour of
the system.
CHAPTER 1. INTRODUCTION 3

Figure 1.2: The geometry of a 2-D collapsible channel (from paper by Luo [48]).

To overcome the inadequacy of the fluid-membrane model on flow in collapsible chan-

nels, Cai and Luo [17] developed a fluid-beam model, which employs a plane strained
elastic beam with large deflection and incrementally linear extension. This model repre-
sents a more realistic and general description of the problem and can be reduced to several
simpler models including the fluid-membrane model under special parameter ranges. Both
asymptotic and numerical approaches are used to study the problem.
Heil investigated the steady deformations of the fully coupled three-dimensional system
in another series of studies [38], [36], [39]. The wall of the tube was modelled as a circular
cylindrical shell and geometrically nonlinear shell theory was used to describe the large non-
axisymmetric post-buckling deformation. The fluid flow was modelled by using lubrication
theory. But some of the assumptions used in the simplification of the fluid equations were
violated, which caused the wall slope at the downstream end to tend to be quite large,
after the buckling, although this model provides a very accurate description of the tube’s
deformation. Heil [37] developed the entirely three-dimensional self-consistent model of the
viscous flow in a collapsible tube by abandoning the small-slope assumption and replacing
the lubrication theory by a solution of the steady three-dimensional Stokes equations which
describe the flow in arbitrary geometries at zero Reynolds number.
In order to develop a more general three-dimensional model requires one ultimately to
extend the Heil model [37] to replace the geometrically nonlinear shell theory by a theory
which can accommodate the large displacements (geometric nonlinearity) and as well as
large strains (material nonlinearity) and also to include the coupling of unsteady, three-
dimensional, nonlinear Navier-Stokes equations for oscillations to arise. Because of the
complexities and large computational requirements for the full three-dimensional solution
of the above problems, such a work, however, is still a daunting task. Although some
attempts have been made, self-excited oscillations are still not yet captured [71], [70], [65].
CHAPTER 1. INTRODUCTION 4

To reduce the complexity we simplify the actual physiological problem by replacing the
transmural (internal minus external) pressure due to the flow in and out the tube by the
external hydrostatic pressure P only. This simplification allows us to avoid tackling the
fully coupled fluid-structure interactions and to focus on investigating the material and
geometrical nonlinearity of the tubes. We consider the simplified problem in 3 stages
(1) Bifurcation analysis of thick-walled circular cylindrical elastic tubes under axial
loading and external pressure.
(2) Nonlinear axisymmetric deformations of an elastic tube under external pressure.
(3) Nonlinear three-dimensional deformations of an elastic tube under external pres-
sure.

We give a brief review of the stability analysis here. To gain a full understanding of
this particular area, we refer to the recent review article by Fu and Ogden [26], in which
they summarized the progress of development of the nonlinear stability analysis of thick
elastic bodies subjected to finite elastic deformations. The stability of elastic shells has
been analyzed over the course of the past century since the initial work on the topic by von
Mises [72], who derived an equation for the buckling pressure of a thin-walled elastic tube.
This gives the pressure as proportional to the cube of the ratio of wall thickness to mean
diameter. Since then buckling of circular cylindrical tubes under external pressure based
has been studied extensively, for instance by Batdorf [7], Nash [54] and Flügge [25]. In
these studies, a simple one-term deflection function was used and the problem was solved
under special boundary conditions. More accurate solutions were obtained by Ho [40],
Sobel [66] and Yamaki [75] for a variety of loading and boundary conditions where the
pre-buckling state was given in terms of membrane theory. The same problem was then
treated by Yamaki [76] but with pre-buckling effects. His key finding was that the mode
number of the most unstable mode increases as the tube length is decreased, and for
a sufficiently long tube mode 2 bifurcation is the most unstable mode. The length of
the tube at which the transition between the higher mode and mode 2 occurs, however,
depends on the thickness ratio; the thicker the tube the shorter the length for which mode 2
becomes the most unstable mode [77]. Good agreement between these studies and various
experiments [74] has led to the buckling prediction for a cylindrical tube being regarded
as a solved problem (at least for thin shells).
With different emphases, related extensive studies on stabilities of circular cylindrical
CHAPTER 1. INTRODUCTION 5

shells have also been carried out. Some of these, concerning geometrically nonlinear vi-
brations and dynamics of circular cylindrical shells, were reviewed by Amabili [1], with
and without fluid-structure interactions. Other recent advances in post-buckling analysis
of thin-walled structures were reported by Kounadis [46]. With a particular interest in
post-buckling behaviour, Heil and Pedley [39] examined the stability of cylindrical shells
under external pressure using a geometrically non-linear shell theory and confirmed that
the mode number of the most unstable mode increases as the tube length is decreased, as
predicted by Yamaki [77]. Heil and Pedley [39] also found that the bifurcation is not signif-
icantly affected by the presence of a full fluid-solid coupling (as long as the critical loading
is the same), although the subsequent post-buckling behaviour can be very different with
and without the internal flow.
In experimental studies for these kinds of problems the tube wall thickness typically
exceeds that which might be appropriate for thin-shell theories [10, 11, 13]. It is there-
fore reasonable to ask if the bifurcation predictions of the classical theories remain valid.
Bertram [10] studied experimentally the effects of wall thickness on the collapse of tubes
and obtained agreement with the results of [74]. In Bertram’s study, wall thickness ratio
h/R values used were 0.38 and 0.5, where h is the thickness and R is the internal radius.
The thick-walled tube problem was also analyzed by Marzo [52] using the finite element
method, and good agreement with the experiments of [10] and [74] was achieved. However,
in [10] and [52] results were presented only for mode 2 bifurcation and for limited values
of the wall thickness. Therefore, it remains unclear how far the bifurcation predictions
of thin-shell theory can be extended to thick-walled tubes, for which nonlinear elastic
deformations can no longer be neglected.
There is also an extensive literature on plastic buckling of circular tubes. Experimental
and modelling aspects of the compression of steel tubes in the plastic regime have been
reviewed in the recent works by Bardi [5] and [6]. They found that the carbon steel tubes
may buckle into different modes as the increase of the external pressure. Figure 1.3 shows
the plastic buckling of circular tubes under compression with axisymmetric collapse and
non-axisymmetric collapse, with mode 2 and 3. This figure can give the reader directly an
idea of the shape of the tube after axisymmetric or non-axisymmetric buckling. We refer
to these papers for references to the relevant literature.
For problems involving nonlinear elastic deformations, a rigorous bifurcation theory
has been established based on the analysis of infinitesimal deformations superimposed on
CHAPTER 1. INTRODUCTION 6

Figure 1.3: (a) Carbon steel tube that developed axisymmetric concertina folding, (b)
mode 2 folding and (c) mode 3 folding of stainless steel tubes (from paper by Bardi [5]).

a known large deformation [28]. Using this theory, Nowinski and Shahinpoor [56] exam-
ined the stability of an infinitely long circular cylinder of neo-Hookean material under
external pressure assuming a plane strain deformation, and Wang and Ertepinar [73] in-
vestigated the stability of infinitely long cylindrical shells and spherical shells subjected to
both internal and external pressure. On the same basis but for different (incompressible,
isotropic) material models Haughton and Ogden [35] examined in some detail the bifur-
cation behaviour of circular cylindrical tubes of finite length under internal pressure and
axial loading.
Bifurcation from a circular cylindrical configuration of a thick-walled tube subject
to combined axial loading and external pressure was investigated on the basis of the
nonlinear theory of elasticity by Zhu et al. [78]. Our work showed that the wall thickness
and aspect ratio play important roles in the occurrence of the most unstable bifurcation
mode. Different from the results based on thin shell theories, which show that higher
modes should occur for shorter tubes, Zhu et al. [78] showed that mode-2 becomes more
persistent for shorter tubes if a suitable nonlinear model is used. This observation was in
agreement with experimental findings on thick-walled tubes subject to external pressure,
in particular those of [10, 11] and [13].
However, a limitation of this work is that the bifurcation analysis was initiated from
a deformed circular cylindrical configuration of an elastic tube with rather special incre-
mental boundary conditions imposed on the ends of the tube. Thus, the results only apply
for the initial bifurcation behaviour, and might preclude realistic post-buckling behaviour
CHAPTER 1. INTRODUCTION 7

involving large displacements near the ends of the tube.

In many engineering and biomechanical applications cylindrical tubes are subject to
external pressures and as a result undergo large (nonlinear) deformations. In another work
we investigated the behaviour of thick-walled tubes with large deformations including large
strains and large displacements under external pressure by deriving the general differential
equations, free of unnecessary assumptions. This is a challenging work due to the lack of
the good and general numerical methods for solving the final fully nonlinear differential
equations we obtained based on the theory of nonlinear finite deformation and the large
computational requirements for simulation of the full three dimensional problems, so ini-
tially, we assume the deformation is axisymmetric. Then, we will try to move forward to
the full three-dimensional problems.
In early engineering approaches to the analysis of this problem it was typically as-
sumed that the material response is linearly elastic, but this led to predictions which were
inaccurate except for very small deformations. It is well known that for biological ma-
terials deformations of the order 50–100% can occur, and in this case a fully nonlinear
problem formulation is essential. However, fully nonlinear material and geometrical anal-
ysis is challenging due to the difficulty of solving such problems. To facilitate solutions
simplifications are often made, such as the adoption of thin shell theories, which have been
successful for describing thin-walled structures [47, 75, 77]. Some researchers have focused
on geometrically nonlinear problems, with small strains but large displacements, and this
approach has often proved to be adequate. Erbay and Demiray [23] considered the finite
axisymmetric deformation of a circular cylindrical tube of neo-Hookean material by using
an asymptotic expansion method. Their perturbation solution is based on the smallness
of the ratio of thickness to inner radius of the tube. Normal and tangential tractions were
applied on the inner surface of the tube but no boundary conditions were considered at
the ends of the tube. Heil [38] and Marzo [52] performed a numerical simulation of the
post-buckling behaviour of tubes under external pressure.
Propagation of finite amplitude waves in fluid-filled elastic or viscoelastic thin-walled
tubes has been investigated [64], [2], [53], and [23]. However, for thick-walled tubes there
are few results available in the literature due to the difficulties arising from the variation
of field quantities with the radial coordinate. Demiray studied weakly nonlinear waves in a
fluid filled thick-walled elastic tube, first using an artificial estimated pressure dependence
[20] on the axial coordinate, which was later improved upon [21].
CHAPTER 1. INTRODUCTION 8

The ability to predict the bifurcation character of the solutions is also an important
practical problem. Negrón-Marrero [55] studied the bifurcation of the axisymmetric hy-
perelastic cylinders subject to nonlinear mixed boundary conditions and found that the
eigenfunctions can be classified into those that are symmetric about the mid-plane, rep-
resenting either necked or barrelled configurations of the cylinder, and those that break
this symmetry. Finite axisymmetric deformations of thick-walled carbon-black filled rub-
ber tubes were also studied experimentally by Beatty and Dadras [9]. They found that
for aspect ratios less than 5 tubes exhibit radially or axially symmetric bulging modes
of deformation, distinct from the familiar Euler buckling that occurs for longer tubes.
Significantly, they found that the experimentally observed critical compression load is
considerably lower than that predicted on the basis of the linear theory.
In chapter 2, we introduce the theory of nonlinear elasticity, which will be used through-
out the thesis.
In chapter 3, the finite element process and techniques that are used in Chapter 5 are
presented. Particular attention and details are provided to introduce the object-oriented
finite element library libmesh, which will be adopted for solving the nonlinear partial
differential equations in Chapters 5 and 6.
In chapter 4, following the analysis of [35], we consider the bifurcation of incompress-
ible, isotropic thick-walled circular cylindrical tubes of finite length when subject to both
axial loading and external pressure. A new feature of the present work is the combination
of finite deformations of thick-walled tubes of hyperelastic material with external pressure
and axial loading.
For the thinner tubes it is found that under external pressure axisymmetric bifurcation
occurs only for 0 < λz < 1, where λz is the principal stretch in the axial direction of the
finite deformation. Moreover, the trend of the bifurcation curves is very different from
that of a tube under internal pressure. Since externally pressurized tubes are particularly
prone to asymmetric bifurcations, we devote most of our effort to the study of asymmetric
bifurcations. The bifurcation modes are characterized by azimuthal mode number m and
the tube length (which can be taken as a proxy for the axial mode number n). The
bifurcation curves for modes m = 1 to m = 4 are presented, and the effects of wall
thickness and the ratio of tube length to external radius on the buckling pressure are also
examined. For the simpler cases, our results are in agreement with the published results
in [52], [74], [10] and [73], and, in particular, with the von Mises equation [72, 74]. We
CHAPTER 1. INTRODUCTION 9

observe that the von Mises equation can only predict the buckling pressure well for thin
shells. By contrast, the general analysis of bifurcation based on 3D finite deformation
elasticity theory presented herein is valid for both thin and thick shells.
In chapter 5, we formulate the fully nonlinear problem of the large axisymmetric de-
formations of thick-walled cylindrical tubes of finite length made of incompressible hy-
perelastic material subject to zero displacements on the ends of the tube and hydrostatic
pressure on the exterior of the lateral surface. The general governing differential equations
that describe the deformation of the tube are derived, with both geometrical and material
nonlinearity included. The corresponding radially symmetric and linear problems are also
examined for the purpose of comparison. The sets of equations are solved numerically
using the object-oriented C++ finite element package Libmesh. Results for tubes with
different aspect ratios are presented to show how the wall thickness and tube length affect
the nonlinear behaviour. The major findings are that for a short tube with smaller aspect
ratio, the nonlinear deformation is characterized by a corner bulging, which changes all
the stress distributions, especially for the shear stress. For longer tubes, the nonlinear
model exhibits higher modes of deformation while for the corresponding linear model only
mode-2 is present. The agreement between the linear and nonlinear models is only good
for small values of the pressure, corresponding to maximum strains of about 5%.
In chapter 6, without any assumptions on the magnitude of the geometrical deforma-
tion or material nonlinearity we derived the general three dimensional governing equations
for the large deformations of a thick-walled tube composed of incompressible isotropic elas-
tic material in both cylindrical polar and Cartesian coordinates. Generally, it is convenient
that we formulate our equations for a circular cylindrical tube based on cylindrical polar
coordinates. However, we note that the expression of deformation gradient F in cylindrical
polar coordinates is much more complex than the one in Cartesian coordinates and this
complexity can be enlarged in the expression of nominal stress S and even the equilibrium
equations. The form of the final equilibrium equations in cylindrical polar coordinates is
also more complicated, with several redundant terms. Both of the complexities will add
difficulty when the numerical discretization of the equation system and computations are
carried out. In order to avoid the complexities in formulation, we prefer to adopt the corre-
sponding Cartesian equation systems, although dealing with the boundary condition may
seem to be not rational compared with an approach based on cylindrical polar coordinates.
The only thing we need do is to get the expression of the unit normal to the internal and
CHAPTER 1. INTRODUCTION 10

external surfaces of the cylindrical tube. The corresponding linear equations in Cartesian
coordinates are also presented for the purpose of comparison with the nonlinear ones.
Results from Chapters 4 and 5 have been published in International Journal of Solids
and Structures [78] and European Journal of Mechanics A solids [79], respectively. Further
results from Chapter 6 are still in preparation and will appear soon.
Chapter 2

Basic equations

In this chapter a brief summary of the static theory of nonlinear elasticity is given, includ-
ing the analysis of deformation, strain, stress and the governing equations of equilibrium,
and a short description of the constitutive equations for a Cauchy elastic material. For
more important details we refer to the relevant literature such as the classic book by Og-
den [60], in which a complete and precise account of the mathematical theory of non-linear
elasticity with application to the analysis of the large elastic deformation of hyperelastic
materials is presented and the book by Fu and Ogden [27], which provides not only fun-
damentals of nonlinear elasticity but also modern topics in this field.

2.1 Deformation

We will deal with deformations of elastic material in which both rotations and stretches
are arbitrarily large, the so-called finite strain theory. In this case, a clear distinction is
necessary to be made between undeformed and deformed configurations of an elastic body.
Consider a deformable continuous body for which we take X to be the position vector of
an arbitrary material point in the reference configuration, denoted by Br . Similarly, in
the current configuration, Bt say, let x be the position vector of the same material point.
Suppose that the deformation from Br to Bt is defined by the vector function χ, if
there is no time dependence we have that (see Fig.2.1) x = χ(X). We assume that χ is
twice-continuously differentiable with respect to position here.
The displacement vector u is defined by

x = X + u. (2.1)

11
CHAPTER 2. BASIC EQUATIONS 12

dA F da
u

X x

O o
(1) Reference configuration (2) Current configuration

Figure 2.1: Reference and current configurations.

Then the deformation gradient tensor F is defined by

F = Gradx. (2.2)

With respect to Cartesian basis vectors Ei in reference configuration and ei in current

configuration, we have
∂(xi ei ) ∂xi
F= ⊗ Ej = ei ⊗ Ej i, j = 1, 2, 3. (2.3)
∂Xj ∂Xj
The local ratio of current to reference volume is

J = det F > 0, (2.4)

and for an incompressible material the constraint

J = det F ≡ 1 (2.5)

must be satisfied for every material point X.

For any non-singular second order tensor F, we note that the tensor can be written
uniquely in the form
F = RU = VR, (2.6)

where R is a proper orthogonal tensor, so that

RRT = RT R = I, (2.7)

where I is the identity tensor. The tensors U and V are positive definite and symmetric,
the so-called right and left stretch tensors, respectively. The eigenvalues of U are the
(strictly positive) principal stretches of the deformation, denoted λi , i = 1, 2, 3. Please
note that λi are also the eigenvalues of V. Then by using (2.7) we can easily get

J = det F = det U = det V, (2.8)

CHAPTER 2. BASIC EQUATIONS 13

In 1839, George Green [18] introduced a deformation tensor called the right Cauchy-Green
deformation tensor or Green’s deformation tensor, which is defined as

C = FT F = U2 . (2.9)

Physically, the tensor C gives us a measure of local change in length of an line element
due to deformation.
It is also useful to note that the Nanson’s formula is given by

nda = JF−T NdA, (2.10)

where dA is the area element of material surface in Br and da is the corresponding area
element in Bt ; See Fig.2.1. n and N are the unit outward normals in the current and
reference configurations, respectively. The connection (2.10) can be used to map from
areas in the current configuration to the corresponding areas in the reference configuration
and vice versa.

2.2 Stress theory and equilibrium

Let t denote surface (contact) force, per unit deformed area, which depend continuously
on x and n.

Theorem 2.2.1 Cauchy’s theorem: if t(x, n) is continuous in x, then there exists a

second-order tensor field σ such that

t(x, n) = σ(x)n, (2.11)

where the tensor σ is also called the Cauchy stress tensor and is independent of n.

The Cauchy stress tensor σ is symmetric, i.e. σ T = σ, and satisfies the Eulerian form of
the equilibrium equation, namely

divσ + ρb = ρa, (2.12)

where ρ is the mass density of the material of the body in current configuration Bt and b
is the body forces, measured per unit volume. a is the acceleration.
We can write the surface force on an area element nda in the current configuration as
following by using (2.10) and Cauchy theorem (2.11)

tda = σnda = JσF−T NdA = ST NdA, (2.13)

CHAPTER 2. BASIC EQUATIONS 14

where the relation of nominal stress tensor S and Cauchy stress tensor σ is given by

S = JF−1 σ. (2.14)

The corresponding nominal stress tensor S, also referred to as the engineering stress,
which, in general, is not symmetric, satisfies

FS = ST FT . (2.15)

Please note that ST is the first Piola-Kirchhoff stress tensor.

The Lagrangian alternative to the Eulerian equilibrium equation (2.12) is

DivS + ρr b = ρr a, (2.16)

where the Div is the divergence operator with respect to X and the mass density ρr is
related to ρ by the mass conservation equation

J = ρr /ρ. (2.17)

Then, in the static case, if there are no body forces the local equilibrium equation for
the body has the (Lagrangian) form

DivS = 0, (2.18)

or, in terms of Cauchy stress,

divσ = 0. (2.19)

2.3 Constitutive law for a Cauchy elastic material

In solid mechanics, a constitutive equation approximates the actual response of a material

to external forces. To be more precise it connects the stresses to strains or stretches.
A simple elastic material is one for which the stress at each material point is dependent
solely on the current state of deformation with respect to an arbitrary reference config-
uration. However, the work done by the stress does in general depend on the path of
deformation and the stress cannot be derived from a scalar potential function.
Neglecting the effect of temperature the constitutive equation for a homogeneous elastic
material can be written as

σ = g(F), (2.20)
CHAPTER 2. BASIC EQUATIONS 15

where g is called response function of the material relative to Br , which is a symmetric

tensor-valued function. We can see that the Cauchy stress σ at an arbitrary material point
X is determined only by the deformation gradient F at this point and doesn’t depend on
the history of deformation.
The principle of objectivity requires that material properties should be independent of
superposed rigid-body motions. This means the constitutive law g must satisfy

g(QF) = Qg(F)QT , (2.21)

for each F and any rotation Q, which is a proper orthogonal second-tensor.

If for all proper orthogonal second-order tensors Q, we have

g(FQ) = g(F), (2.22)

then the material is said to be isotropic relative to Br . In essence, this means the material
properties have no preferred direction.
In equation (2.22), with Q replaced by RT and use of polar composition (2.6), we get

σ = g(F) = g(VRRT ) = g(V). (2.23)

Using material objectivity (2.21) combined with the definition of isotropy (2.22) and
(2.23), we obtain

g(QFPT ) = Qg(FPT )QT = Qg(F)QT = Qg(V)QT , (2.24)

then choose P = QR, and we have

g(QVQT ) = Qg(V)QT , (2.25)

which then shows that g is an isotropic tensor function of V. It can be shown that the
Cauchy stress σ may be written in the form
3
X
σ= σi v(i) ⊗ v(i) , (2.26)
i=1

where

σi = φ0 + φ1 λi + φ2 λ2i i = 1, 2, 3. (2.27)

φi = φi (I1 , I2 , I3 ) and the invariants are defined by

1
I1 = tr(C), I2 = [I12 − tr(C2 )], I3 = det C. (2.28)
2
CHAPTER 2. BASIC EQUATIONS 16

2.4 Green elastic material

A Green elastic or hyperelastic material is an ideal special case of a Cauchy elastic material
for which a strain-energy function exists. The observed material behaviour of rubber, filled
elastomers and biological tissues are often described by the hyperelastic idealization. The
constitutive relation of such a material can be defined as isotropic, incompressible, non-
linearly elastic and generally independent of strain rate.
The strain-energy function W (F) is given by
∂
W (F) = Jtr(σL), (2.29)
∂t
where the velocity gradient tensor L, is defined as

L = gradv, (2.30)

where v is the velocity vector. Physically, W (F) is a measure of the work per unit refer-
ence volume done by stress as a result of deformation and is independent of the path of
deformation.
We also have
µ ¶
∂ ∂W
W (F) = tr Ḟ , (2.31)
∂t ∂F
combined with Ḟ = LF, we can get
µ ¶ µ ¶
∂ ∂W ∂W
W (F) = tr LF = tr F L . (2.32)
∂t ∂F ∂F
Comparison of this with (2.29) shows that stress tensor σ can be written in terms of W (F)
as
∂W
σ = J −1 F , (2.33)
∂F
or in component form (Cartesian coordinates),
∂W
σij = J −1 Fik . (2.34)
∂Fjk
Note that the components of ∂W/∂F are defined by the convention
µ ¶
∂W ∂W
= . (2.35)
∂F ij ∂Fji

Using the connection (2.14) between the nominal stress S and the Cauchy stress tensor σ,
it follows that
∂W
S= , (2.36)
∂F
CHAPTER 2. BASIC EQUATIONS 17

and the component form is

∂W
Sij = . (2.37)
∂Fji
Definition 2.4.1 Since W is a scalar function objectivity requires that it is unaffected by
a superimposed rigid-body rotation after deformation, i.e.

W (QF) = W (F) (2.38)

for all rotations Q and for each deformation gradient F.

Definition 2.4.2 For a hyperelastic material which is isotopic relative to Br , W is un-

affected by rotations in Br (prior to deformation), such that

W (FPT ) = W (F) (2.39)

for all rotations P.

Using the definitions of objectivity and isotropy of W (F), we could deduce

W (QVQT ) = W (V), (2.40)

for all rotations Q. This means that W is an isotropic scalar function of V. Thus, the strain
energy function W has all the properties associated with the isotropic scalar function, i.e.
it is expressible as a function of the principal invariants I1 , I2 , I3 or equivalently, as a
symmetric function of the principal stretches λ1 , λ2 , λ3 . So, we have

W (λ1 , λ2 , λ3 ) = W (λ1 , λ3 , λ2 ) = W (λ3 , λ1 , λ2 ). (2.41)

2.4.1 Stress-deformation relations in terms of invariants

We regard W as a function of the invariants I1 , I2 , I3 , defined in equation (2.28), i.e.

W (I1 , I2 , I3 ). Then, we could express nominal stress S as
3
X ∂W ∂Ii
S= . (2.42)
∂Ii ∂F
i=1

Using the connection of (2.14) between the nominal stress S and the Cauchy stress
tensor σ, we could easily get the corresponding Cauchy stress
3
X ∂W ∂Ii
σ= J −1 F . (2.43)
∂Ii ∂F
i=1

where the derivatives are

∂I1 ∂I2 ∂I3
= 2FT , = 2I1 FT − 2FT FFT , = 2I3 F−1 . (2.44)
∂F ∂F ∂F
CHAPTER 2. BASIC EQUATIONS 18

2.5 Internal constraints

To simplify the constitutive response models and also represent a good first approximation
to the actual material behaviour, some form of local constraints are used, such as incom-
pressibility or inextensibility. In mathematical theory, the question is how the constraints
influence the evaluation of the stress tensor.
Suppose the deformation is constrained by the single scalar function

C(F) = 0. (2.45)

Differentiation of (2.45) with respect to time gives

µ ¶
∂C
Ċ ≡ tr Ḟ = 0. (2.46)
∂F

Compared with the stress power defined by (2.31), we could accommodate the constraint
in the stress-deformation relation by adding an arbitrary scalar multiple of C(F), without
affecting the stress power i.e.

W (F) + qC(F), (2.47)

where the hydrostatic pressure q functions is a Lagrange multiplier, in general, q is inde-

pendent of F and dependent on X.
The Cauchy stress tensor σ and nominal stress tensor S defined by (2.33) and (2.36)
respectively are modified to

∂W ∂C
Jσ = F + qF , (2.48)
∂F ∂F

and

∂W ∂C
S= +q . (2.49)
∂F ∂F

For a material constrained by incompressibility, we have

C(F) = J − 1 = λ1 λ2 λ3 − 1 = 0. (2.50)

To ensure the incompressibility of an elastic material, we can replace strain energy function
W (F) by

W (F) − p(det F − 1). (2.51)

CHAPTER 2. BASIC EQUATIONS 19

Thus, the nominal stress S and Cauchy stress σ are given by

∂W
σ=F − pI, (2.52)
∂F

and

∂W
S= − pF−1 . (2.53)
∂F

2.6 Example strain-energy functions for isotropic elastic ma-

terial

The neo-Hookean material model is given by

1
W (λ1 , λ2 , λ3 ) = µ(λ21 + λ22 + λ23 − 3), (2.54)
2

where µ(> 0) is a material constant referred to as the shear modulus of the material. This
is an extension of Hooke’s law for the case of large deformations and can be applied to
plastics and rubber-like substances. However, the neo-Hookean material model usually
provides sufficient accuracy for materials under moderate straining up to 30-70%.
The Mooney-Rivlin material model, a generalization of the neo-Hookean model, is
defined by

1 1
W (λ1 , λ2 , λ3 ) = µ1 (λ21 + λ22 + λ23 − 3) − µ2 (λ21 λ22 + λ22 λ23 + λ23 λ21 − 3), (2.55)
2 2

where µ1 (≥ 0) and µ1 (≤ 0) are material constants such that µ1 − µ2 = µ(> 0).

For complex materials such as rubbers, polymers, and biological tissue subject to even
larger deformation, more sophisticated models are necessary. The Ogden material model,
which was developed by Ray W. Ogden in 1972, like other hyperelastic material models,
assumes that the material behaviour can be described by a strain energy density function,
from which the stress-deformation relationships can be derived.
The Ogden material model is given by
N
X µn αn
W (λ1 , λ2 , λ3 ) = (λ + λα2 n + λα3 n − 3) (2.56)
αn 1
n=1

where µn and αn are material constants and satisfy the constraint as follows
N
X
µn αn = 2µ, n = 1, 2, 3, ..., N, (2.57)
n=1
CHAPTER 2. BASIC EQUATIONS 20

where N is a positive number. Note that choosing appropriate material constants µn ,

αn and N , the Ogden model can be reduced to Mooney-Rivlin or neo-Hookean material
model.
For more details of strain-energy functions in terms of principle stretches we refer to
Ogden [57], [59].

2.7 Incremental deformations

We now consider a deformation relative to a given reference configuration Br , defined by

x = χ(X), (2.58)

then superimpose an infinitesimal deformation on the known deformation χ(X), such that

0 0
x = χ (X). (2.59)

The infinitesimal displacement of the body due to this change is denoted by

0
δx = x − x. (2.60)

Here we assume the displacement δx is small enough so that the terms of order |δx|2 and
higher order can be neglected in comparison with those of order |δx|. δx is referred to as
an incremental deformation from the configuration described by χ(X).
The corresponding change of deformation gradient due to the incremental displacement
δx is then given by

Gradδx = δGradx = δF. (2.61)

Using the definition of the differentiation of a scalar function of a tensor, described in

Section 4.2.8 in the book by Ogden [60], and the Taylor series, we obtain the change of a
scalar function as
µ ¶ µµ 2 ¶ ¶
∂φ 1 ∂ φ
δφ(F) = tr δF + tr δF δF + higher orders. (2.62)
∂F 2 ∂F2
In terms of Cartesian components the first term in equation (2.62) is defined by
µ ¶ µ ¶
∂φ ∂φ
tr δF = δFji , (2.63)
∂F ∂Fji
and the second term in component form is defined by
µµ 2 ¶ ¶
1 ∂ φ 1 ∂2φ
tr δF δF = δFlk δFji . (2.64)
2 ∂F2 2 ∂Fji ∂Flk
CHAPTER 2. BASIC EQUATIONS 21

Let φ = J and neglect the terms of second and higher order in the definition of (2.62),
then δJ, the change in the determinant of deformation gradient can be written as

δJ = Jtr(F−1 δF) (2.65)

The increment of a tensor function G can be treated in a similar manner, so that

µ ¶
∂G 1 ∂2G
δG = δF + δF δF + higher orders, (2.66)
∂F 2 ∂F2
with Cartesian components
µ ¶
∂G ∂Gij
δF = δFkl , (2.67)
∂F ij ∂Fkl
and
µµ ¶ ¶
∂2G ∂ 2 Gij
δF δF = δFmn δFkl . (2.68)
∂F2 ij ∂Fkl ∂Fmn
Then the linear approximation of the nominal stress increment may be written in
accordance with (2.66)

δS = AδF (2.69)

where the notation A, referred to as the tensor of first-order elastic moduli associated
with the conjugate pair (S, F), is defined by
∂S
A= , (2.70)
∂F
which is also called the tensor of fixed-reference moduli.
If we take the reference configuration to coincide with the current configuration Bt ,
then the increment of the deformation gradient relative to current configuration takes the
form

δF0 = δFF−1 , (2.71)

where the subscript zero indicates a quantity evaluated in Bt .

We also have the corresponding connection

δS0 = J −1 FδS, (2.72)

Then the resulting elastic moduli are referred to as the instantaneous moduli and the
relation between instantaneous moduli and fixed-reference moduli is given by (for details,
see [60])

A0ijkl = J −1 Fiα Fkβ Aαjβl . (2.73)

Please note that we also introduce the alternative notation ẋ = δx = η, Ḟ = δF, Ṡ = δS

and the instantaneous moduli A0 = B in Chapter 4.
Chapter 3

Finite element method, libmesh

library

3.1 Finite element nonlinear analysis in solid mechanics

3.1.1 Introduction

The finite element method, which originated from the need for finding approximate solu-
tions to complicated structural analysis and elasticity problems in civil and aeronautical
engineering, is now an important and indispensable tool in scientific research, engineering
analysis and design, such as in the thermal, electromagnetic, solid, fluid, and structural
working environments.
An alternative way of solving partial differential equations is the finite difference
method. Compared to the finite difference method, the most attractive feature of the
finite element method is its ability to handle complicated geometries with relative ease.
While finite difference methods can be very easy to implement, in general, the accuracy
of a finite element method approximation is often more precise than in the corresponding
finite difference method.
To alleviate difficulties in solving problems with localized features that are not effi-
ciently resolved by mesh refinement, the extended finite element method, also known as
the generalized finite element method was developed in 1999 by Belytschko and collabo-
rators. This method has been used to model the propagation of various discontinuities,
such as cracks and material interfaces.
Over the past 20 years, meshfree methods have been developed to facilitate simulations

22
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 23

in problems where the ability to handle discontinuities and singularities, large deforma-
tions, advanced materials is needed. For example, the melting of a solid or the freezing
process can be simulated using meshfree methods.

3.1.2 Procedure for finite element solution with libmesh

Discussions on the finite element method in detail are obviously beyond the scope of this
work. However, we will summarize the particular techniques and the processes related to
the finite element analysis in Chapter 5.
For discretization of the differential equations, the weighted residual-Galerkin method
has been used. This method is frequently adopted in the finite element literature since it
usually (but by no means always) leads to symmetric matrices [80]. In essence the original
shape functions are used as weighting functions for the approximations used in the integral
formulations (see Appendix 1).
The displacement-based finite element procedures are not sufficiently effective for the
analysis of incompressible materials, and mixed finite element models have therefore been
adopted to obtain good solution accuracy [8], as illustrated in Fig. 3.1. This figure
shows that, for a 6-node triangle element, the displacements are interpolated using the six
nodes and the pressure, which is a Lagrange multiplier coming from the incompressibility
constraint, see equation (2.5), is interpolated by using only three corner nodes. On the
other hand, for 9-node quadrilateral elements, the displacements are interpolated using
nine nodes and the pressure using 4 corner nodes only. Mathematically, the linear and
bilinear pressure interpolations are used respectively for the above two cases, i.e.

p = p0 + p1 x + p2 y,
p = p0 + p1 x + p2 y + p3 xy. (3.1)

In general, once a mathematical model has been developed, the numerical finite element
procedure can be followed to solve the governing equations approximately. In the following
we summarize the implementation of these procedures with an object-oriented parallel
finite element library, which will be covered in the next section. The procedure is depicted
in Fig. 3.2. Discretization of the governing equations leads to expressions of the element
b e using
stiffness matrices Ke . We could obtain the modified element stiffness matrices K
boundary conditions. Instead of forming all elements first and then assembling them,
we will construct elements one at a time in a loop, and immediately merge them into
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 24

(a) (b)
Node with displacement variables

Node with displacement and pressure variables

Figure 3.1: Mixed elements: (a) 6-node Triangle (b) 9-node Quadrilateral.

Finite element solution Mathematical model

Discretized equations Differential equations
Element type/Mesh Geometry
Global stiffness Matrix Kinematics
Linear/nonlinear solver Material law
Stress recovery Boundary conditions

Element
Ke Impose l
K
e
Assembler
Stiffness Boundary
Matrices Conditions
l R
K
Post Nodal Equation
processing displacement solver

Figure 3.2: Procedures in a Finite element program.

b The residual vector R can be obtained based on K.

the Global matrix K. b Then, the

linear/nonlinear solver can be called to solve the final algebra equations. Once we get
the values of the displacements and the pressure, post processing can be carried out to
obtain the stresses or stretches. We could use libmesh to assemble the system and call
the equation solvers to find the approximate solution. However, since libmesh is not a
black box tool and still underdeveloped now, we have to create our own mesh files in 2D
problems or use Tetgen (an open-source tetrahedral mesh generator) to mesh the tube in
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 25

our 3D problem. Also, we have to write code for the purpose of post processing.

3.1.3 Stress recovery

It is very important to evaluate the stresses at each node for all the elements in engineering
application. A quite straightforward approach is to evaluate the stresses at the nodes of a
given element by substituting the natural coordinates of the nodes to the shape functions
then using the connection between stress and displacements. Another approach is to
evaluate the stresses at Gauss integration points and then extrapolate to the element
nodes [22]. The latter approach provides more accurate stress values for quadrilateral
elements, since the best accuracy for gradients and stresses is obtained at the Gauss
points. We will adopt the second approach for the calculation of stresses in Chapter 5.
Other variables such as the principal stresses/stretches are also evaluated in a similar
way; see details in Chapter 5. The detailed explanation of this method is given as below.
However, for triangle elements both of the above approaches give similar results. Note
that this discussion is mainly based on the course materials from the Department of
Aerospace Engineering Sciences, University of Colorado at Boulder; for details, see the
website: http://www.colorado.edu/engineering/cas/courses.d/IFEM.d/.
Taking a 4-node quadrilateral element for example, the four Gauss points, denoted as
0 0 0 0
1 , 2 , 3 , 4 are composed of the “Gauss element”, as illustrated in Fig.3.3. The natural
coordinates of both the nodes of the quadrilateral element and the nodes of Gauss element
0
are listed in Table 3.1. The “Gauss element”, denoted by (e ), is also a 4-node quadrilateral,
0 0
with its coordinates denoted by ξ and η . The connections between two sets of coordinates
are
0 √ 0 √
ξ = ξ / 3, η = η / 3. (3.2)

An arbitrary scalar quantity u can be approximated by

4
X 0
0 0 (e ) 0 0 0
u(ξ , η ) = Ni (ξ , η )ui , (3.3)
i=1
0
(e )
where Ni , i = 1, 2, 3, 4 are the shape functions, defined by
(e )
0
1 0 0
N1 = (1 − ξ )(1 − η ),
4
0
(e ) 1 0 0
N2 = (1 + ξ )(1 − η ),
4
0
(e ) 1 0 0
N3 = (1 + ξ )(1 + η ),
4
0
(e ) 1 0 0
N4 = (1 − ξ )(1 + η ). (3.4)
4
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 26

When the stresses at Gauss element points are evaluated, the extrapolation procedure can
0 0
be done. For example, to extrapolate u to node 1, we replace its coordinates (ξ , η ) =
√ √
(− 3, − 3) into equation (3.3).

Table 3.1: Natural coordinates of Quadrilateral nodes.

0 0 0 0
Corner nodes ξ η ξ η Gauss nodes ξ η ξ η
√ √ 0 √ √
1 −1 −1 − 3 − 3 1 −1/ 3 −1/ 3 −1 −1
√ √ 0 √ √
2 +1 −1 + 3 − 3 2 +1/ 3 −1/ 3 +1 −1
√ √ 0 √ √
3 +1 +1 + 3 + 3 3 +1/ 3 +1/ 3 +1 +1
√ √ 0 √ √
4 −1 +1 − 3 + 3 4 −1/ 3 +1/ 3 −1 +1

0
Figure 3.3: (a) 4-node quadrilateral element (e); (b) Gauss element (e )

Figure 3.4 shows Gauss elements for 8-node and 9-node quadrilaterals and a 6-node
triangle. Extrapolation in these higher order elements can be evaluated in a similar way
easily. However, it is a process demanding great caution for the implementation of this
stress recovery technique compatibly with libmesh.
In finite element analysis, it is an assumption that the elements must be complete
and compatible. The compatibility condition requires the displacements and their mth
derivatives are continuous across the adjacent element for a C m variational problem [8]. In
the analysis of a plate bending problem, for example, the transverse displacement u is the
only unknown variable. The transverse displacement u and its derivatives ∂u/∂x, ∂u/∂y
are continuous. This continuity condition can’t guarantee the continuity of the stresses
calculated at the same node of adjacent elements. This indicates some necessary process
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 27

Figure 3.4: Gauss elements for high order quadrilaterals and triangles (a) 9-node element
with 3 × 3 Gauss rule (b) 9-node element with 3 × 3 Gauss rule (c) 6-node element with
3-interior rule.

of stress averaging is needed to smooth and improve the accuracy of the stresses. The
stress averaging might be followed in two ways in practice:
1. Unweighted averaging: assign the same weight to all elements that share a common
node.
2. Weighted averaging: the weight assigned to element contributions depends on the
stress component and the element geometry and possibly the element type.
For the problems in Chapter 5 the unweighted averaging is chosen to compute the
averaged nodal stresses σij as well as other nodal variables such as the principal stretches
λi .

3.2 Libmesh library

The open-source finite element library libmesh [44] provides a software framework for
parallel adaptive finite element simulations of partial differential equations using arbitrary
unstructured discretization. It is also integrated with third party software packages such as
(1) PETSc and LASPack for the solution of linear systems on both serial and parallel plat-
forms, (2) METIS and ParMETIS for mesh partitioning, (3) Triangle and Tetgen for mesh
generation. libmesh has been developed at The University of Texas at Austin in the CFD
Lab since March 2002. The libmesh library itself is not tied to any particular application.
The simulations presented in Chapter 5 were performed by the author using application
codes built on top of this framework. The libmesh library is coded by the object-oriented
C++ programming language. In the following section, an overview of Object-Oriented
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 28

scientific computing is presented, which is a key to gaining better understanding of the

libmesh framework.

3.2.1 Object-Oriented Scientific Computing

The C++ programming language provides many useful features for simulating complex sci-
entific computing problems [15], [16], [68], such as abstraction, encapsulation, inheritance,
polymorphism. A critical feature missing in Fortran 90 is the template techniques, which
allows C++ programmers to build portable, reusable code and to dramatically improve
the efficiency of the evaluation of complex expressions involving user-defined data types.
In recent years the performance of the C++ programming language has improved. Many
high performance, object-oriented scientific softwares have been developed using C++.
There are two distinct paradigms for implementing software algorithms, namely procedure-
oriented and object-oriented approaches.
The procedural approach has dominated numeric computation and scientific comput-
ing for decades. One of the most popular procedural programming language is Fortran. In
this approach a sequence of algorithmic steps operates on some set of data structures to
implement a given algorithm. In consequence the data storage and procedure implemen-
tation are intimately related. Suppose a standard array were used to store the individual
elements of a vector, for example. If for some reason, a linked-list would be a more efficient
data structure, due to dynamic insertion and removal of elements, for example, then it
would require substantial changes to all codes which use such a vector.
On the other hand, object-oriented approaches provide user-defined classes which define
the attributes and the behaviours of a particular data type. The class concept is a tool
that can be used to create new data types. Within a given class, data and function
members can be declared as either public, protected, or private in order to explicitly enforce
encapsulation. A significant feature of classes is encapsulation. As a result, the actual data
is separated from operations which are performed on the data. Considering the vector
example again, object-oriented programming allows that the specific data structure used
to store the elements of the vector can be completely encapsulated within an object and the
codes which use such an object don’t need to have any access to this data structure. Then,
if the algorithmic implementation or data storage techniques of an object are changed for
some reason, the codes using such an object need not to change. For this and many
other reasons, object-oriented programming has been used widely since the mid-1990s
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 29

to build more maintainable and extensible software. The application of object-oriented

techniques in scientific numerical simulations has been slow but much effort is being made
to implement high performance object-oriented scientific software [3], [4], [44].

3.2.2 libmesh: Adaptive mesh refinement scheme

Figure 3.5: Element refinement hierarchy for a 2D quadrilateral mesh (from PhD thesis
by Kirk [45]).

One of the main features of libmesh is the support for adaptive mesh refinement. In an
adaptive refinement procedure, when a solution on a given mesh is obtained, the error of
this solution will be estimated to get local error indicators which can be used as the criteria
for selective local mesh refinement. Two main categories of procedures for the adaptive
refinement of the finite element solutions are the h-refinement and p-refinement. In h-
refinement the elements are changed in size, as illustrated in Fig. 3.5; some of the elements
are made larger and others are smaller. By contrast, p-refinement keeps the size of the
elements but increases the order of the polynomial used in their definition [19], [80]. Both
of the refinement approaches are provided in the libmesh library. In a refinement process
used by libmesh a new set of children elements is created from the parent elements through
a linear map, which is provide by an “embedding matrix”. On the other hand, in the
coarsening process all the children of a given element are removed and the parent element
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 30

is re-activated. For details of the error indicators and refinement criteria, we refer to the
work done by Kelly [43]. Figure 3.6 shows the solution to the nonlinear model (described
in Chapter 5) for an adapted quadrilateral mesh using the KellyErrorEstimator class in
libmesh.

Figure 3.6: Adaptive mesh refinement on a rectangle domain.

3.2.3 Data structure of libmesh

Most of this section is based on the libmesh web page: http://libmesh.sourceforge.net/,

the paper [44] and the dissertation by Kirk [45].
The Mesh class provides a description of a geometric entity. A mesh is composed of
elements and nodes, which are stored in the mesh. These data are encapsulated by abstract
classes which provide an interface for a variety of possible implementations. To access the
particular subset or all the nodes and elements in the mesh, the user just needs to create a
node/element iterator object. In addition, this class provides functions for implementing
mesh input/output in various formats, including GMV format from Los Alamos National
Labs, TetGen, Tecplot, Exodus II from Sandia National Labs, and GMSH.
The abstract base class Elem provides an interface for implementation of a geometric
element. The derived classes, such as Hex8 support the actual operation and calculations
for a given geometric elements via virtual functions. Figure 3.7 shows a simplified inher-
itance diagram, in which a Cell is an abstract Elem in three-dimensions and a cell could
be a tetrahedron, a pyramid, a prism, or a hexahedron. The concrete subclass Hex8 is
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 31

an element composed of 8 nodes in three dimensions. The user can conveniently obtain
the number of the nodes of all geometric element types by calling the virtual function
n_nodes() on an Elem pointer. These virtual functions allow for defining a new element
type by the user without affecting the external application programming interface; for
example, the original code used to return the number of nodes of a given element type.

Figure 3.7: A simplified inheritance diagram for the Dof Object class (from PhD thesis by
Stogner [67]).

In a classic finite element data structure, the element connectivity is usually given in
terms of the node indices, while in the libmesh library the Elem class stores pointers to
the nodes to which the element is connected. This approach can enable the element to
determine the location of its nodes with a single pointer dereference. Elements also contain
pointers to their face neighbours and their parent or child elements. When at least one
side of the element is on the physical boundary of the domain, it means the element has a
Null neighbour with a Null pointer added into the array of the pointers to its neighbour.
It is convenient to apply boundary conditions by finding all the elements with a Null
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 32

pointer neighbour.
The abstract System class contains information related to a set of differential equations
that might be simulated. Several concrete systems are derived, such as LinearImplicitSystem,
NonlinearImplicitSystem which will be used to solve the linear and nonlinear sets of
equations in Chapters 5 and 6. Note that a system is uniquely tied to a particular mesh;
in a simulation multiple meshes are used, and then multiple systems are needed.
The base class NonlinearSolver provides a uniform interface for nonlinear solvers in
packages like PETSc. PETSc, a portable and extensible library for scientific computing,
is the underlying parallel linear solver used in this work, which was developed in the
Mathematics and Computer Science Division at Argonne National Laboratory [4].
Although libmesh offers all the standard geometric element types, such as triangles,
quadrilaterals, tetrahedra, hexahedra, prisms and pyramids, it can’t automatically gener-
ate meshes for complex geometries but is only limited to simple geometries like a rectangle,
a circle and a cube. For the two dimensional problem described in more details in Chapter
5, a symmetric mesh is necessary to obtain the correct solution, since the geometry of
the tube section and the boundary conditions, including pressure on lateral surface of the
tube and zero-displacement end conditions, are all symmetric. We have written several
mesh files (for details, see Appendix) in XDA format to store the coarse symmetric mesh
data and then using uniform refinement loops to obtain the final mesh. For the three
dimensional problems in Chapter 6, we use the mesh generator TetGen [31] to generate
tetrahedral meshes, which then can be read into libmesh. The C++ code used to obtain
the tetrahedral mesh (as illustrated in Fig. 3.8) is provided in Appendix 2.
CHAPTER 3. FINITE ELEMENT METHOD, LIBMESH LIBRARY 33

X Y

0.5

z
0
-1 -1

-0.5 -0.5

0 0
y x
0.5 0.5

1 1

Figure 3.8: A 3D tube discretized by tetrahedral element.

Chapter 4

Bifurcation

In this chapter we restrict our attention to buckling of circular cylindrical tubes under
external pressure.
In Section 1 we summarize the necessary constitutive equations that describe finite
elastic deformations, while in Section 2 these are specialized to the circular cylindrical
geometry of a thick-walled tube that maintains its circular cylindrical shape under axial
extension and external pressure. The equations that describe a general (three-dimensional)
incremental deformation superimposed on the deformed circular cylindrical tube are then
given in Section 3. The three coupled partial differential equations governing the in-
cremental displacement components are highlighted in Section 4 along with the relevant
incremental boundary conditions. Based on an appropriate Ansatz for the displacement
components the equations reduce to coupled ordinary differential equations, for the solu-
tion of which a numerical scheme is then described. In Section 5 the numerical method
is used in respect of a specific material model in order to obtain details of the onset of
bifurcation in either an axisymmetric or asymmetric mode.

4.1 The elastic constitutive law and strain-energy function

We consider the material body to be composed of an elastic material, whose properties

are described in terms of a strain-energy function, which we denote by W = W (F) per
unit reference volume. Here we confine attention to incompressible materials, so that the
stress-deformation relation is given by either

∂W
S= − pF−1 , (4.1)
∂F

34
CHAPTER 4. BIFURCATION 35

where p (an arbitrary hydrostatic stress) is a Lagrange multiplier associated with the
constraint (2.5), or
∂W
σ=F − pI, (4.2)
∂F
where I is the identity tensor.
Here we take the material to be isotropic, so that W depends on F only through
the principal stretches λi , i = 1, 2, 3, and is a symmetric function of the stretches. We
therefore represent W in the form W = W (λ1 , λ2 , λ3 ), and, for an incompressible material,
the constraint (2.5) may be written in terms of the stretches as

λ1 λ2 λ3 = 1. (4.3)

Moreover, (4.2) can be decomposed on principal axes as

∂W
σi = λi − p, i = 1, 2, 3 (no summation), (4.4)
∂λi

σi , i = 1, 2, 3, being the principal Cauchy stresses.

For subsequent convenience it is useful to regard W as a function of just two indepen-
dent stretches, λ1 and λ2 say, and to introduce the notation Ŵ defined by

Ŵ (λ1 , λ2 ) = W (λ1 , λ2 , λ−1 −1

1 λ2 ). (4.5)

It then follows from (4.4) that the principal stress differences can be written

∂ Ŵ ∂ Ŵ
σ1 − σ3 = λ1 , σ2 − σ3 = λ2 . (4.6)
∂λ1 ∂λ2

4.2 The circular cylindrical configuration

We now consider a thick-walled circular cylindrical tube with reference geometry described
by
A ≤ R ≤ B, 0 ≤ Θ ≤ 2π, 0 ≤ Z ≤ L, (4.7)

where R, Θ, Z are cylindrical polar coordinates, A and B are the inner and outer radii,
respectively, and L is the length of the tube. This is depicted in Fig. 4.1(a).
The initial deformed configuration of the tube, under the action of axial loading and
external pressure, is assumed also to be circular cylindrical, with geometry described by

a ≤ r ≤ b, 0 ≤ θ ≤ 2π, 0 ≤ z ≤ l, (4.8)
CHAPTER 4. BIFURCATION 36

A
R
B

Z
Ĭ O

(a)

a
b r
a
o
ș L

(b)

Figure 4.1: The circular cylindrical tube in its reference configuration (a) and deformed
configuration when subject to axial load and external pressure (b).

where r, θ, z are cylindrical polar coordinates, a and b are the internal and external radii,
respectively, and l is the length. Since the material is incompressible, the deformation is
described by the equations

r2 = a2 + λ−1 2 2
z (R − A ), θ = Θ, z = λz Z, (4.9)

where λz is the axial extension ratio (or axial stretch), which is uniform.
We use e1 , e2 , e3 to denote the unit basis vectors corresponding to the coordinates
θ, z, r, respectively. For the considered deformation, since the material is isotropic, these
CHAPTER 4. BIFURCATION 37

define the principal directions of both the stretch tensor U and the Cauchy stress σ.
Let λ1 , λ2 , λ3 denote the corresponding principal stretches and σ1 , σ2 , σ3 the associated
principal Cauchy stresses, which are given by (4.4). From the incompressibility constraint
together with (4.9), we have

r
λ2 = λz , λ1 = ≡ λ, λ3 = (λ1 λz )−1 . (4.10)
R

For the symmetric configuration considered here, the only equilibrium equation not
satisfied trivially is
dσ3
r + σ3 − σ1 = 0, (4.11)
dr
and we have the associated boundary conditions

 0 on r = a
σ3 = (4.12)

−P on r = b.

Using Ŵ , as defined in (4.5), (4.6)1 , and the definitions (4.10), integration of (4.11)
and application of the boundary conditions (4.12) yields
Z b
dr
P =− λŴλ . (4.13)
a r

On application of the connections r = λR and (4.9) this may be re-written with λ as the
integration variable in the form
Z λb
Ŵλ
P = dλ, (4.14)
λa (λ2 λz − 1)

where
a b
λa = , λb = . (4.15)
A B
We note here that if there is, additionally, an internal pressure, Pi > 0 say, then the
left-hand sides of (4.13) and (4.14) are replaced by P − Pi . Thus, the effect of an internal
pressure can be captured by taking P < 0 in the above formulas, this corresponding to a
radial external tension on r = b.

4.3 Incremental equations

Detailed derivation of the incremental equations can be found in [35] for a thick-walled and
[34] for a thin-walled tube (see [32] and [33] for corresponding results for spherical shells).
Here we provide a summary of the main results needed for our analysis. A superposed
CHAPTER 4. BIFURCATION 38

dot signifies an increment in the quantity concerned, and a subscript 0 indicates that the
quantity to which it is attached is calculated with respect to the deformed configuration
as reference configuration. First, let ẋ(X) denote the incremental displacement vector,
and then define u(x) through u(x) = u(χ(X)) = ẋ(X). Note that u was used earlier for
the displacement in equation (2.1) in Chapter 1, which does not appear in this chapter so
there is no conflict of notation. Next, introduce the notation η defined by

η = Ḟ0 ≡ ḞF−1 = gradu. (4.16)

The incremental form of the incompressibility condition can then be written

trη = 0. (4.17)

The increment of the constitutive law (4.1) has the form

Ṡ = AḞ − ṗF−1 + pF−1 ḞF−1 , (4.18)

where A is the elasticity tensor with components defined by

∂2W
Aαiβj = . (4.19)
∂Fiα ∂Fjβ
When the reference configuration is updated to the current configuration this becomes

Ṡ0 = Bη + pη − ṗI, (4.20)

where I is again the identity tensor and B is the 4th-order tensor of instantaneous elastic
moduli, whose (Cartesian) components are related to those of A by

Bpiqj = Fpα Fqβ Aαiβj . (4.21)

For an incompressible isotropic elastic material the non-vanishing components of B

referred to the principal axes of σ can be written (see, for example, [58])

Biijj = Bjjii = λi λj Wij , (4.22)

λi Wi − λj Wj 2
Bijij = λi , λi 6= λj , (4.23)
λ2i − λ2j
Bijji = Bjiij = Bijij − λi Wi , i 6= j, (4.24)

where Wi = ∂W/∂λi , Wij = ∂ 2 W/∂λi ∂λj .

The incremental form of the equilibrium equation (2.18) is Div Ṡ = 0 and when updated
it becomes
div Ṡ0 = 0, (4.25)
CHAPTER 4. BIFURCATION 39

the incremental counterpart of (2.19).

For the problem to be considered in the following sections we shall be making use of
the pressure boundary condition, which, referred to the original reference configuration,
may be written
ST N = −P F−T N, (4.26)

where N is the unit outward normal vector to the boundary of the body in the reference
configuration and P is the pressure on the boundary per unit area of the deformed con-
figuration. On taking the increment of (4.26) and updating to the deformed configuration
we obtain
ṠT T
0 n = P η n − Ṗ n, (4.27)

which is the form of incremental boundary condition that we shall use.

We now specialize (4.25) to circular cylindrical coordinates based of the underly-
ing solution discussed in Section 3. The curvilinear coordinates are ordered so that
(x1 , x2 , x3 ) = (θ, z, r). Then, we have, in component form,

Ṡ0ji,j + Ṡ0ji ek · ej,k + Ṡ0kj ei · ej,k = 0, i = 1, 2, 3, (4.28)

with summation over indices j and k from 1 to 3, where the subscript j (= 1, 2, 3) following
a comma represents the derivatives (∂/r∂θ, ∂/∂z, ∂/∂r). The only non-zero components
of ei · ej,k are
1 1
e1 · e3,1 = , e3 · e1,1 = − . (4.29)
r r
Referred to the cylindrical polar axes the incremental displacement u is written in
terms of its components (v, w, u) as

u = ve1 + we2 + ue3 . (4.30)

Then, from the definition η = gradu we obtain the component matrix of η referred to the
axes in question as  
(u + vθ )/r vz vr
 
 
[η] =  wθ /r wz wr  , (4.31)
 
(uθ − v)/r uz ur
where the square brackets indicate the matrix of components of the enclosed quantity and
the subscripts (r, θ, z) signify standard partial derivatives.
The incompressibility condition (4.17) can now be given explicitly as

trη ≡ ur + (u + vθ )/r + wz = 0. (4.32)

CHAPTER 4. BIFURCATION 40

4.4 Asymmetric bifurcations and numerical methods

We now substitute (4.20), (4.31), (4.32) and the expressions for the components of Bijkl
into (4.28) to obtain

0
ṗθ = (rB3131 + B3131 )(uθ + rvr − v)/r + (B1111 − B1122 − B2112 )(uθ + vθθ )/r
+ B2121 rvzz + B3131 rvrr + (B1133 − B1122 − B2112 + B3113 )urθ , (4.33)

0
ṗz = (rB3232 + B3232 )(uz + wr )/r + B1212 (wθθ − ruz )/r2 + B3232 wrr
+ (B2222 − B1221 − B1122 )wzz + (B2233 + B3223 − B1221 − B1122 )urz , (4.34)

0 0
ṗr = (rB1133 − rB2233 − B1111 + B1122 )(u + vθ )/r2 + B1313 (uθθ − vθ )/r2 + B3223 wrz
+ (B1331 + B1133 − B2233 )vrθ /r + (B3333 − B2233 )urr + B2323 uzz
0
+ (rB3333 + rp0 − rB2233
0
+ B3333 − 2B2233 + B1122 )ur /r. (4.35)

On the cylindrical boundaries we apply the specialization of (4.27) to the present

situation, with the inner boundary free of incremental traction and the outer boundary
subject to pressure P . Taking Ṗ = 0 in (4.27) we then have, for i = 1, 2, 3,

0 on r = a
Ṡ03i = (4.36)

P η3i on r = b.

At the ends of the tube we apply the incremental boundary conditions

u = v = 0, Ṡ022 = 0 on z = 0, l. (4.37)

This means that the ends of the tube are constrained so that no incremental rotation or
radial displacement is allowed, while the axial component of traction is of dead-load type.
To solve the equations, we assume that the solution takes the form

u = f (r) cos mθ sin αz, v = g(r) sin mθ sin αz, 

(4.38)

w = h(r) cos mθ cos αz, ṗ = k(r) cos mθ sin αz, 

where m = 0, 1, 2, 3, . . . is the azimuthal mode number, m = 0 corresponding to an

axisymmetric solution. Substitution into the incompressibility condition (4.32) then yields

rf 0 (r) + f (r) + mg(r) − αrh(r) = 0. (4.39)

Also, on inserting (4.38) into (4.33)–(4.35) and using (4.39) to eliminate h(r), we obtain
CHAPTER 4. BIFURCATION 41

three coupled equations for f (r), g(r) and k(r), namely

0
(rB3131 + B3131 + B1111 − B1122 − B2112 )mf (r) + (B1133 − B1122 + B3113 − B2112 )mrf 0 (r)
0
+ [rB3131 + B3131 + m2 (B1111 − B1122 − B2112 ) + α2 r2 B2121 ]g(r)
0
− (rB3131 + B3131 )rg 0 (r) − B3131 r2 g 00 (r) − mrk(r) = 0, (4.40)

0
[rB3232 − B3232 + m2 B1212 − α2 r2 (rB3232
0
+ B3232 − B1212 + B1122 + B1221 − B2222 )]f (r)
0
− [rB3232 − B3232 − m2 B1212 − α2 r2 (B2222 − B2233 − B3223 )]rf 0 (r)
0
− (rB3232 + 2B3232 )r2 f 00 (r) − B3232 r3 f 000 (r)
0
+ [rB3232 − B3232 + m2 B1212 + α2 r2 (B2222 − B1122 − B1221 )]mg(r)
0
− (rB3232 − B3232 )mrg 0 (r) − B3232 mr2 g 00 (r) + α2 r3 k(r) = 0, (4.41)

0 0
(rB1133 − rB2233 − B1111 + B1122 + B3223 − m2 B1313 − α2 r2 B2323 )f (r)
0
+ (rB3333 + rp0 − rB2233
0
+ B3333 − 2B2233 + B1122 − B3223 )rf 0 (r)
+ (B3333 − B2233 − B3223 )r2 f 00 (r)
0 0
+ (rB1133 − rB2233 − B1111 + B1122 + B3223 − B1313 )mg(r)
+ (B1133 − B2233 + B1331 − B3223 )mrg 0 (r) − r2 k 0 (r) = 0. (4.42)

Next, on substituting the expression for u from (4.38) in the boundary condition
(4.37)1 , we deduce that
α = nπ/(λ2 L), (4.43)

where n = 1, 2, 3, . . . is the axial mode number. The boundary conditions for v are then
automatically satisfied. It is therefore clear that the behaviour for different mode numbers
n can be captured, equivalently, by varying the length L. Thus, in what follows it suffices
to set n = 1 and to consider L as a parameter that reflects either changes in the axial
mode number or changes in length.
¿From equations (4.40)–(4.42), we can express f 000 (r), g 00 (r) and k 0 (r) in terms of f (r),
f 0 (r), f 00 (r), g(r), g 0 (r) and k(r), and hence we write the equations as a first-order system
in the compact form
dy
= G(y, r), (4.44)
dr
where y = (y1 , y2 , y3 , y4 , y5 , y6 )T , G = (G1 , G2 , G3 , G4 , G5 , G6 )T ,

y1 = f (r), y2 = f 0 (r), y3 = f 00 (r), y4 = g(r), y5 = g 0 (r), y6 = k(r), (4.45)

and

G 1 = y2 , G2 = y3 , G4 = y5 , (4.46)
CHAPTER 4. BIFURCATION 42

while G3 , G5 , G6 are lengthy expressions obtained by rearranging equations (4.40)–(4.42)

and are not listed here.
In the same notation, the components of the incremental pressure boundary condition
(4.27) are given as

my1 + y4 − ry5 = 0, 




(α2 r2 + m2 − 1)y1 + ry2 + r2 y3 = 0, (4.47)





(B1133 − B2233 )(y1 + my4 ) + (B3333 − B2233 + λ3 W3 )ry2 − ry6 = 0,
each of which must hold on both r = a and r = b. To obtain these use has been made of
the conditions σ3 = λ3 W3 − p = 0 on r = a and σ3 = λ3 W3 − p = −P on r = b, and we
have set Ṗ = 0 on r = b.
To solve the system of first-order ordinary differential equations (with three indepen-
dent solutions), we choose starting values at r = a for three independent solutions given
by    
y11 (a) y12 (a) y13 (a) 100
   
 1   
 y4 (a) y42 (a) y43 (a)  =  0 1 0  , (4.48)
   
y61 (a) y62 (a) y63 (a) 001

where, for each entry yij (a) in (4.48), subscripts i = 1, 4, 6, correspond to dependent
variables in (4.44) while the superscript j refers to the jth set of initial values (j = 1, 2, 3).
Substituting each set of the initial values, that is each column of the matrix (4.48),
into the boundary conditions (4.47) for r = a, we obtain
   
y21 (a) y22 (a) y23 (a) a11 my21 (a) a13
   
 1   
 y3 (a) y32 (a) y33 (a)  =  a21 −my21 (a)/a −y23 (a)/a  , (4.49)
   
y51 (a) y52 (a) y53 (a) m/a 1/a 0
where, for conciseness, we have introduced the notations
B2233 − B1133 1
a11 = , a13 = ,
a(B3333 − B2233 + λ3 W3 ) B3333 − B2233 + λ3 W3
1 − m2 − a2 α2 − ay21 (a)
a21 = ,
a2
all terms being evaluated for r = a.
Equations (4.48) and (4.49) together give the initial values for equations (4.44). This
initial value problem is solved numerically using the Adams-Moulton method (Gerald and
Wheatley, 1984), with Predictor and Corrector given by
h
Predictor: yn+1 = yn + (55Gn − 59Gn−1 + 37Gn−2 − 9Gn−3 ), (4.50)
24
CHAPTER 4. BIFURCATION 43

h
Corrector: yn+1 = yn + (9Gn+1 + 19Gn − 5Gn−1 + Gn−2 ), (4.51)
24
where h = (b − a)/ω is the step size and ω is the iteration number. Note that the Adams-
Moulton method requires four sets of initial values at previous steps. These are calculated
using the fourth-order Runge-Kutta method. Each method has local errors of O(h5 ).
The solutions can be expressed as a linear combination of the three independent solutions
y1 , y2 , y3 . Thus,
y = C1 y1 + C2 y2 + C3 y3 , (4.52)

where yi = (y1i , y2i , y3i , y4i , y5i , y6i )T , i = 1, 2, 3.

Bifurcation may occur if there exist constants C1 , C2 , C3 , at least one of which is non-
zero. For purposes of numerical computation in Section 6 we shall specialize to a particular
strain-energy function, for which B1133 = B2233 = 0. On introducing this specialization
and substituting (4.52) into the boundary conditions (4.47), we obtain three equations for
C1 , C2 , C3 , namely

[my1i (b) + y4i (b) − by5i (b)]Ci = 0, 




[(m2 + α2 b2 − 1)y1i (b) + b(y2i (b) + by3i (b))]Ci = 0, (4.53)




[b(B3333 + λ3 W3 )y2i (b) − by6i (b)]Ci = 0, 

evaluated for r = b, in each of which there is summation over the index i from 1 to 3. Thus,
the bifurcation criterion is obtained by the vanishing of the determinant of coefficients of
C1 , C2 , C3 , viz.
¯ ¯
¯ 1 1 1 2 2 2 3 3 3 ¯
¯ my1 (b) + y4 (b) − by5 (b) my1 (b) + y4 (b) − by5 (b) my1 (b) + y4 (b) − by5 (b) ¯
¯ ¯
¯ 1 1 2 1 2 2 2 2 3 3 2 3 ¯
¯ M y1 (b) + by2 (b) + b y3 (b) M y1 (b) + by2 (b) + b y3 (b) M y1 (b) + by2 (b) + b y3 (b) ¯ = 0,
¯ ¯
¯ ¯
¯ 1 1
bN y2 (b) − by6 (b) 2 2
bN y2 (b) − by6 (b) 3 3
bN y2 (b) − by6 (b) ¯
(4.54)
again with all terms evaluated for r = b, where M = m2 + α2 b2 − 1 and N = B3333 + λ3 W3 .
Substituting the equation

b2 = a2 + λ−1 2 2
z (B − A ), (4.55)

i.e. equation (4.9)1 with R = B, into (4.54), we obtain an equation for the value of a that
satisfies the bifurcation criterion (4.54). The corresponding bifurcation pressure can be
obtained from (4.14).
CHAPTER 4. BIFURCATION 44

4.5 Numerical results and discussion

In the experiments of [74] and [10] silicone rubber tubes were used, and the numerical
results of [52] were compared with experimental data for two thick-walled collapsible tubes
reported by [10]. It is therefore appropriate to employ a strain-energy function that has
been used extensively for fitting data on experiments for a wide range of rubberlike solids.
Specifically, we apply the foregoing theory to the strain-energy function given by
3
X
W = µr (λα1 r + λα2 r + λα3 r − 3)/αr , (4.56)
r=1

where µr and αr , r = 1, 2, 3, are material constants (see, for example, [60]). Using the
incompressibility condition (4.3) and the energy function Ŵ (λ1 , λ2 ) defined by (4.5), we
have
3
X
Ŵ (λ1 , λ2 ) = µr (λα1 r + λα2 r + (λ1 λ2 )−αr − 3)/αr . (4.57)
r=1
For the numerical calculations we use the material constants given by

α1 = 1.3, α2 = 5.0, α3 = −2.0,

µ∗1 = 1.491, µ∗2 = 0.003, µ∗3 = −0.023, (4.58)

as in [33], where µ∗r = µr /µ, r = 1, 2, 3, and µ is the shear modulus of the material in the
reference configuration given by (see, for example, [57])
3
X
2µ = µr αr . (4.59)
r=1

Representative values of the aspect ratios of the tube are taken as L/B = 1, 2.5, 5, 10,
and for numerical purposes, without loss of generality, we set B = 1 and change the value
of the inner radius A to vary the thickness of the tube. Two thickness ratios are considered,
namely, A/B = 0.85 (thinner tube) and A/B = 0.5 (thicker tube).
The qualitative nature of the results presented below are not unduly sensitive to the
choice of material parameters in (4.56), and there are also many other forms of strain-
energy function that could equally well be used to produce similar qualitative behaviour.

4.5.1 Equilibrium pressure curves

The dependence of the non-dimensional pressure P ∗ = P/µ on the circumferential stretch

λa is illustrated in Fig. 4.2(a) in respect of the strain-energy function (4.57) with material
constants (4.58) and for A/B = 0.85 and several values of λz . Figure 4.2(a) shows that
CHAPTER 4. BIFURCATION 45

initially the external pressure increases slowly in order to compress the tube radially as λa
is reduced from 1. Thereafter, there is a plateau where a significant increase in pressure
does not produce significant further radial deformation of the tube. This trend becomes
more pronounced as the value of λz increases. This graph should be compared with the
pressure-area (internal cross sectional area of the tube) diagram, also known as the “tube
law” and most commonly used for collapsible tubes [24]. Although the tube law is based
on the post-buckling behaviour of tubes it doesn’t take account of axial forces and bending
moments.
The equal pressure curves corresponding to P ∗ = 0, 0.5, 1 are plotted in (λz , λa ) space
for A/B = 0.5 and 0.85 in Fig. 4.2(b), again using equation (4.14), except for P ∗ = 0, for
which we have the connection
λ2a λz = 1, (4.60)

which is independent of the wall thickness ratio A/B. We observe that at least for the
range of values of λz and λa considered, the equal pressure curves for the thicker tube
(A/B = 0.5) lie above those for the thinner one (A/B = 0.85), indicating that to obtain
the same deformation more pressure is required for the thicker tube, as should be expected.

4.5.2 Axisymmetric bifurcation

First, we consider axisymmetric modes of bifurcation, corresponding to m = 0 in (4.38).

We set the longitudinal mode number n to be 1 and in Fig. 4.3 we plot axisymmetric
bifurcation curves for L/B = 2.5, 5, 10 and 20 and A/B = 0.85. In this case, as well
as curves for an external pressure, curves for an internal pressure are shown in order
to compare with the results of [35]. With reference to the remarks on internal pressure
following equation (4.14), we recall that the effect of internal pressure is captured by
taking P ∗ < 0 here. It can then be seen that for a tube subjected to internal pressure our
results coincide with those in [35] except for a factor 2, which means the curves in [35] for
L/B = 2x are the same for those here with L/B = x.1
When the tube is under external pressure (P ∗ > 0), we note that the axisymmetric
bifurcation curves all intersect the curve P ∗ = 0 in the region 0 < λz < 1, which means that
axisymmetric bifurcation cannot occur for tubes with A/B = 0.85 subjected to external
pressure and axial extension (i.e. when λz > 1). In other words, under external pressure,
1
Private communication with Dr. Haughton confirms that there is a factor of 2 missing in eq. (61)
of [35].
CHAPTER 4. BIFURCATION 46

P*
5

1 5 4 3 2 1

Λa
0.2 0.4 0.6 0.8 1

(a)

Λa

1.4 1 P* =0
2 P* =0.5 AB=0.5

1.2 3 P* =0.5 AB=0.85

4 P* =1 AB=0.5
1.0 5 P* =1 AB=0.85

0.8

0.6 1

0.4 2

4 3
0.2 5

0.0 Λz
0 1 2 3 4 5 6

(b)

Figure 4.2: Plot of (a) the dimensionless pressure P ∗ = P/µ against λa for A/B = 0.85
and λz = 1, 2, 3, 4, 5, and (b) equal pressure curves in (λz , λa ) space for P ∗ = 0, 0.5, 1, with
A/B = 0.85 (dashed curves) and A/B = 0.5 (continuous curves).

axisymmetric bifurcation only occurs when a tube is axially compressed. This is not the
case for tubes under internal pressure [35].

4.5.3 Asymmetric bifurcation

Since for tubes under external pressure, axisymmetric bifurcations do not occur when the
tube is extended, we focus on asymmetric bifurcations henceforth.
CHAPTER 4. BIFURCATION 47

Λa
4

3.5

3 20
2.5
2.5 10 P* <0
5
2

1.5

1 P* =0

0.5
P* >0
Λz
1 2 3 4 5

Figure 4.3: Plots of the axisymmetric bifurcation curves for mode n = 1 with aspect ratios
L/B = 2.5, 5, 10, 20 and A/B = 0.85. The dashed curve corresponds to the zero pressure
curve P ∗ = 0.

Thinner tube

In this section, all results are for the thinner tube A/B = 0.85. ¿From equation (4.43),
we recall that either axial mode number n or length of the tube L can be varied to obtain
equivalent results. We therefore set n = 1 and choose different values of the length L,
and only azimuthal modes corresponding to m = 1, 2, 3, 4 are considered. Therefore, in
the following, the mode number referred to is always the azimuthal mode number m. We
restrict attention to m 6 4 because higher mode number bifurcations are not usually
observed in collapsible tube experiments. In any case, we have found that higher modes
produce results very similar to those for m = 4. The asymmetric bifurcation curves
are plotted using the bifurcation criterion (4.54) and the numerical method discussed in
Section 5.
Figure 4.4 shows the mode 1 asymmetric bifurcation curves for L/B = 1, 2.5, 5, 10 and
both internal and external pressure. For P ∗ < 0 (tubes under internal pressure), the results
here are again in agreement with those of [35], with the factor 2 difference indicated earlier,
and we do not discuss this case further. For P ∗ > 0 (tubes under external pressure), we see
that as the axial stretch λz is increased towards 1, along the equal pressure curve P ∗ = 0
the value of λa at bifurcation decreases as the value of L/B increases from 2.5 to 10. This
confirms the intuitive expectation that longer tubes buckle more easily than shorter ones.
CHAPTER 4. BIFURCATION 48

Λa
4

2.5 5
3.5 1 10

2.5

2 P* <0

1.5

1
P* =0
P* >0
0.5

Λz
0.5 1 1.5 2 2.5

Figure 4.4: Mode m = 1 asymmetric bifurcation curves for L/B = 1, 2.5, 5, 10 and A/B =
0.85 in (λz , λa ) space. The dashed curve is the equal pressure curve P ∗ = 0.

Λa
1.6

1.4

1.2
1 LB=1
2 LB=2.5
1
3 LB=5
4 LB=10
0.8
P* =0

0.6

4
0.4 3
2
1
0.2

Λz
1 2 3 4 5 6

Figure 4.5: As in Fig. 4.4 but for azimuthal mode number m = 2.

In the region of axial extension, the tube with L/B = 1 bifurcates slightly more readily
into mode 1 than the longer tubes. Figure 4.4 also shows that the tube can bifurcate into
mode 1 for small axial compression (values of λz less than, but close to, 1). The value of
λa at bifurcation seems to increase rapidly for λz below 1 (i.e. when the tube is axially
compressed). However, under axial extension (λz > 1), bifurcation into mode 1 requires
a relatively larger pressure than in axial compression and the corresponding value of λa
becomes very small, as does the internal radius of the tube.
The mode 2 asymmetric bifurcation curves are shown in Fig. 4.5. It is interesting
CHAPTER 4. BIFURCATION 49

Λa
1.6

1.4
1 LB=1
2 LB=2.5
1.2
3 LB=5
4 LB=10
1
P* =0

0.8

0.6

4
0.4 3
2
1
0.2

Λz
1 2 3 4 5 6

(a)

Λa
1.6
1 LB=1
1.4 2 LB=2.5
3 LB=5
1.2 4 LB=10
P* =0
1

0.8

0.6

4
0.4 3
2
1
0.2

Λz
1 2 3 4 5 6

(b)

Figure 4.6: As in Fig. 4.4 but for (a) m = 3 and (b) m = 4.

to see that the bifurcation pressure for longer tubes (L/B > 5) approaches zero. Thus,
although the bifurcation pressures required in the region of axial compression are similar
for mode 1 and mode 2, much less pressure is required to achieve the mode 2 bifurcation
in the region of axial extension. Figure 4.5 also shows that the mode 2 bifurcation does
not depend significantly on the length of the tube unless the tube is very short (with L/B
about 1).
Similar bifurcation behaviour is found for modes m = 3 and m = 4, as illustrated in
Fig. 4.6. Compared with mode 2, the mode 3 and mode 4 curves are closer to (further
from) the equal pressure line P ∗ = 0 for tubes with L/B = 1 (L/B = 10), and hence
CHAPTER 4. BIFURCATION 50

Λa
1.2
1 m=1

2 m=2
1.0
3 m=3

4 m=4
0.8
P* =0

0.6

0.4 4
3
2

0.2

0.0 Λz
0 1 2 3 4 5 6

(a)
Λa
1.2
1 m=1

2 m=2
1.0
3 m=3

4 m=4
0.8
P* =0

0.6

2
0.4 3
4

0.2

0.0 Λz
0 1 2 3 4 5 6

(b)

Figure 4.7: Asymmetric bifurcation curves for m = 1, 2, 3, 4 and A/B = 0.85 in (λz , λa )
space: (a) L/B = 1; (b) L/B = 5.

the shorter tubes become more sensitive to a change in the external pressure for higher
mode numbers, while for longer tubes, mode 2 become the most unstable mode. The
differences in these modes can be seen more clearly in Fig. 4.7. Note that compared with
higher modes, the mode 1 curve is much further from the P ∗ = 0 curve, especially as axial
extension is increased. This means that unless the tube is slightly compressed, a much
greater pressure is required for a tube to buckle into mode 1 than into higher modes. This
trend is even stronger for the longer tubes.
CHAPTER 4. BIFURCATION 51

Thicker tube

To illustrate the influence of different mode numbers on the behaviour of thicker tubes,
we plot the bifurcation curves for m = 1, 2, 3, 4 in Fig. 4.8 for A/B = 0.5 separately for
each value L/B = 1 and L/B = 5. In Fig. 4.8(a), for L/B = 1, it can be seen that the
bifurcation behaviour for the thicker tube is similar to that for thinner tube, i.e. curves of
modes 2, 3, 4 are closer to each other than that for mode 1. Thus, under extension the tube
may bifurcate into any of the modes 2, 3, 4 but a relatively larger pressure is needed for
mode 1 to be activated. Two major differences are observed between thinner and thicker
tubes. One is that the mode 2, 3, 4 curves are more separated for the thicker tube, the
other is that for axial compression (λz < 1) the lower modes occur first, while for axial
extension, mode 2 becomes the preferred mode for all values of λz . This is consistent with
experimental observations and classical thin shell theory but is not so obvious for thinner
tubes.
Figure 4.8(b) shows corresponding results for L/B = 5. The curves for modes 2, 3, 4
do not intersect. Compared with the L/B = 1 tube, the separations of the curves for
m = 2, 3, 4 are relatively large. The mode 1 curve has one point of intersection with each
of the other higher mode curves. In the region of axial extension, as the external pressure
increases, bifurcation occurs first in mode 2, followed by modes 3, 4 and 1 successively.
For modes 3 and 4, the bifurcation values of λa (larger than 1) along the equal pressure
curve P ∗ = 0 for L/B = 5 are larger than those for L/B = 1.

Bifurcation pressure

Since mode 2 is the most widely observed mode in tube collapse experiments (Bertram,
1987), we show the mode 2 bifurcation pressure against L/B in Fig. 4.9 for both A/B = 0.5
and A/B = 0.85 for comparison, with λz = 1 in each case. It can be seen that the
curves tend to flatten when L/B > 4. This suggests that, for longer tubes, wall thickness
rather than tube length is more important in determining the magnitude of the bifurcation
pressure. As a result, the value of the bifurcation pressure P ∗ for A/B = 0.85 is much
smaller than that for A/B = 0.5, and this will be discussed further later in this section.
It should be noted that different vertical scales are used for the two plots.
To see the change of the bifurcation pressure with wall thickness and to compare our
results with those in the literature (Bertram, 1987; Marzo et al, 2005; Weissman and
Mockros, 1967) we use the reference wall thickness H = B − A and the parameters D, Q
CHAPTER 4. BIFURCATION 52

Λa Λa

2.0 2.0

4
1.5 3 1.5
2 2
1

1
1.0 1.0

P* =0 P* =0

0.5 0.5

0.0 Λz 0.0 Λz
0 1 2 3 4 5 0 1 2 3 4 5

(a) (b)

Figure 4.8: Asymmetric bifurcation curves for m = 1, 2, 3, 4 and A/B = 0.5 in (λz , λa )
space: (a) L/B = 1; (b) L/B = 5.

1.4 0.07

A/B=0.5
A/B=0.85 0.06
1.2

0.05
1

P* A/B=0.85
P* A/B=0.5

0.04
0.8
0.03

0.6
0.02

0.4
0.01

0.2 0
1 2 3 4 5
L/B

Figure 4.9: Plot of P ∗ = P/µ at bifurcation (mode m = 2) against L/B for A/B = 0.5
(continuous curve, left-hand scale) and A/B = 0.85 (dash-dot curve, right-hand scale) and
λz = 1.

and Pk , defined by
2(B − A) EH 3
D= , Q= , (4.61)
ln(B/A) 12(1 − ν 2 )
and µ ¶3
Q 2E H
Pk = 3
= , (4.62)
(D/2) 3(1 − ν 2 ) D
CHAPTER 4. BIFURCATION 53

P
Pk
4

10

3
34

H
0
0.0 0.2 0.4 0.6 0.8 D

Figure 4.10: Mode 2 bifurcation pressure plotted in dimensionless form as P/Pk against
H/D for L/B = 10 (dashed curve) and L/B = 34 (continuous curve) and λz = 1.005.

where, in the context of classical elasticity, E is Young’s modulus and ν is Poisson’s ratio.
Here, D denotes the logarithmic mean diameter and Q is the flexural rigidity of the tube
wall. The pressure P is non-dimensionalized by dividing by Pk .
Using the bifurcation criterion (4.54) combined with equations (4.14), (4.61)1 and
(4.62), we obtain the mode 2 bifurcation pressure shown in Fig. 4.10, plotted with P/Pk
against H/D. We see that the thinner tube begins to bifurcate at a pressure close to
the theoretical value P/Pk = 3 in the thickness range of 0.05 ∼ 0.4 in agreement with
von Mises’ prediction obtained from classical linear elasticity thin shell theory (von Mises,
1914). We emphasize again that our results are obtained from the incremental equations
based on the full 3D theory of nonlinear elasticity, which provide the exact linearized
bifurcation theory of elasticity, and our calculations are valid for underlying finite elastic
deformations. To compare with Bertram’s experimental results [10] and the numerical
results of [52], the parameter L/B = 34 was used here. In fact, our results indicate that,
for tubes with L/B = 10 and L/B = 34, when 0.05 < H/D < 0.4 the values of P/Pk
are in the range 2.9 ∼ 3.2. This explains why von Mises’ prediction is confirmed by
many different experiments and numerical simulations (Bertram, 1987; Marzo et al., 2005;
Weissman and Mockros, 1967). In [10] and [52], only some limited values of P/Pk were
presented for a set of given values of H/D. Likewise, in [74], results were only presented
for 0 < H/D < 0.25. Here, the bifurcation pressure is shown for a much wider range of
H/D. It is interesting to note that the bifurcation pressure does not change significantly
CHAPTER 4. BIFURCATION 54

for tubes with thickness ratio 0.05 < H/D < 0.4.
However, our results show that towards the two ends of the H/D axis, the values of
the bifurcation pressure for mode 2 differ from the classical prediction. For H/D < 0.05,
the values of P/Pk are larger than 3. The shorter the tube, the greater the increase. For
L/B = 34, P/Pk = 3.24 at H/D = 0.01 and for L/B = 10, it increases to 11.5 (see Fig.
4.10 and the B/H = 50 curve in Fig. 4.12). This discrepancy may be because in classical
thin shell theory (Yamaki, 1984) the prebuckling state was assumed to be a membrane
stress state. When H/D < 0.05 and L/B < 34, neglect of the curvature of the deflected
surface caused by external pressure can lead to serious error (von Mises, 1914). However,
von Mises’ formula Pcollapse = 3Pk is sufficiently accurate for shells with L/B > 34 (see
page 73 in Yamaki, 1984). For H/D > 0.4, the curves for L/B = 10 and L/B = 34 almost
coincide. The bifurcation pressure P/Pk drops below 3 as H/D increases, and decreases
to 1 when H/D = 0.8. Caution is required with the physical interpretation of this result,
since Pk is cubic H/D, which increases much faster than P as H/D is increased from 0.4.
This trend can also be seen clearly in Fig. 4.13(a). In physical terms, a greater bifurcation
pressure is still required to buckle the thicker tube, as expected, even though the ratio
P/Pk is smaller.

Very short tubes

To illustrate further the dependence on tube length we now investigate briefly bifurcation
of very short cylinders under axial compression and tension. Figure 4.11(a) presents bi-
furcation curves in (λz , λa ) space for tubes with L/B = 0.5 and A/B = 0.5. Transition
from low to high mode occurs in the range of axial compression at an intersection point
where λz ≈ 0.62. When λz < 0.62, modes 1, 2, 3 occur first, while for λz > 0.62, the mode
m = 4 becomes the most unstable one. Referring back to Fig. 4.8(b) for L/B = 5 we
see that, by contrast, there is no intersection point among curves for m = 2, 3, 4 and the
mode 2 curve is above the others in the whole range of λz except in the short interval
0.90 < λz < 0.95 where the mode 2 curve is below that for mode 1. Axial extension does
not affect the order of the bifurcation modes for either of the tubes with L/B = 1 and
L/B = 5. The parameter L/B therefore plays a major role in the transition from high
to low modes, which is also found for tubes with A/B = 0.85. The results represented
in Fig. 4.11(a) are converted into the plots of P/Pk against λz in Fig. 4.11(b) by use of
(4.61) and (4.62). Figure 4.11(b) shows that the P/Pk curve for mode 1 increases rapidly
CHAPTER 4. BIFURCATION 55

Λa P
Pk

25 1
2.0

20 2

1.5

15

1.0 3

P* =0 10

0.5
5
4
3
1 2
0.0 Λz 0 Λz
0 1 2 3 4 5 6 0 1 2 3 4 5

(a) (b)

Figure 4.11: Asymmetric bifurcation curves for m = 1, 2, 3, 4, L/B = 0.5 and A/B = 0.5
(a): in (λz , λa ) space; (b) in (λz , P/Pk ) space;

and monotonically, while for each mode 2, 3, 4 there is a pressure maximum, occurring at
λz = 1.05, 0.90, 0.80, respectively. Tubes subjected to sufficiently large axial compression
or tension tend to bifurcate easily, while for 0.8 < λz < 1.05 bifurcation requires a larger
pressure. We can therefore conclude that either a large axial compression or axial tension
reduces the axial stiffness of the cylinders.

The most unstable mode

To find the most unstable modes for different lengths and wall thicknesses, similarly to the
predictions of classical thin shell theory (Yamaki, 1984), we plot the critical bifurcation
curves in Fig. 4.12. It is seen that for a thin shell, B/H = 50, the results are in excellent
quantitative agreement with those of Yamaki (1984, figure 2.12, for boundary condition
S4). There exists only a small discrepancy due to the slightly different boundary conditions
used here. In other words, if the wall is thin, then higher modes are more unstable for
shorter tubes. However, as the wall thickness is increased, the critical higher modes
become fewer, and mode 2 becomes more and more dominant. Eventually, for B/H < 2
and L/B > 1.2 it remains the only bifurcation mode. For instance, in the range of
4 < L/B < 10, a thin tube with B/H = 50, bifurcates into the m = 2 mode, whereas
thick-walled tube with B/H = 6.67, bifurcates into the m = 3 mode. In the context of axial
compression of steel cylinders undergoing plastic deformation a very similar distribution
CHAPTER 4. BIFURCATION 56

100
80
60
5

40

4
20 B/H=50
3
P/Pk
3 3

B/H=6.67

B/H=2
2 2
3 2

L/B→ ∞ P/Pk =3.0

B/H=1.58
2 L/B→ ∞ P/Pk =2.43
0 1 2
10 10 10
L/B

Figure 4.12: Bifurcation pressure plotted in dimensionless form as P/Pk against L/B
for B/H = 50 (black curves), B/H = 6.67 (red curves), B/H = 2 (dashed curves),
B/H = 1.58 (blue curves) with different mode numbers and λz = 1.

of bifurcation modes was found by [5] experimentally and [6] analytically. Apart from the
type of material behaviour, this differs from the present analysis since we are considering
external pressure rather than axial compression and we have fixed λz = 1 in Fig. 4.12.
Figure 4.7(a) shows that for tubes with A/B = 0.85 (equivalent to B/H = 6.67) under
external pressure and axial extension, the higher modes are more unstable. Another
interesting phenomenon is that the thicker the tube the smaller the value of L/B at which
the curve flattens. The curves for tubes with B/H = 50, 6.67, 2 show that as L/B → ∞,
P/Pk approaches 3.0, which is in agreement with the thin shell theory prediction. But for
the very thick tube with B/H = 1.58, P/Pk approaches 2.43. The bifurcation pressure for
thick tubes with H/D > 0.4 drops below 3.0 (see also Fig. 4.10).

4.5.4 Discussion

In this chapter, we have investigated the nonlinear buckling behaviour of thick-walled

circular cylinder tubes under external pressure combined with axial loading. Our study is
particularly useful in determining the buckling of thick-walled tubes, which is beyond the
limit of validity of thin shell theory. This work has been conducted with a background in
mind of the bifurcation behaviour of collapsible tubes conveying internal flow. Although
CHAPTER 4. BIFURCATION 57

we note that the essential difference between this study and studies by the collapsible
tube flow community [36, 39] [13] is that no fluid-structure interactions are considered.
Here, the (external) pressure is acting as a (prescribed) static load, which contrasts with
the strong viscous pressure when an internal flow is present. However, in the context of
critical buckling, it has been found that these different mechanisms (static pressure load
or flow-induced pressure load) lead to similar results except that a substantially higher
pressure drop is required to achieve the same level of collapse for the static load case [39].
The most interesting finding is that for wall thickness ratios A/B greater than about
0.5, mode 2 seems to be the dominant critical buckling mode unless the tubes are extremely
short (e.g., L/B . 1.2). This is different from the predictions of classical thin shell theory
[77], but agrees with the fact that in many thick-walled tube experiments, in particular
those of [10, 11] and [13], only mode 2 buckling has been observed regardless of the tube
length used. The fact that in experiments the prevailing mode is mode 2 cannot be
fully explained by thin shell theory. This is because when fluid-structure interaction is
involved, the effect of the fluid flow is to increase the viscous pressure drop, which induces
an additional compressive load at the downstream end of the tube. As a result, only the
compressed downstream part of the tube actually participates in the buckling, which is
then similar to the buckling of a short tube [39]. If the thin shell theory were to be valid,
this would induce the buckling to occur in a higher mode. The reason why this didn’t
happen in the experiments is that, for thicker tubes, mode changes no longer happen, and
long thick tubes were used in experiments [11, 13]. As illustrated in Fig. 4.12, for long
thick tubes, only mode 2 occurs. As indicated above, our study shows that if A/B is
greater than about 0.5, then the critical buckling mode will remain as mode 2 except for
very short tubes.
Although the von Mises formula is derived for thin-walled tubes, experimental mea-
surements have shown that it also predicts the bifurcation pressure for thick-walled tubes
reasonably well [74]. Our results show that this is because the bifurcation pressure P/Pk
is insensitive to the change of wall thickness H/D for the range of 0.05 < H/D < 0.4. If
the tube is sufficiently thin or sufficiently thick, then the von Mises formula is no longer
accurate, and P/Pk actually increases in the thin wall extreme, and decreases in the thicker
wall region.
In order to have a more direct comparison with the Weissman and Mockros experi-
ments, we plot the bifurcation pressure in terms of P against H/D in Fig. 4.13. This
CHAPTER 4. BIFURCATION 58

is obtained using the bifurcation criterion (4.54) combined with equation (4.14) and the
equation
E
µ= , (4.63)
2(1 + ν)
where (for an incompressible material) ν = 0.5. The value E = 300 psi (= 2.07 MPa)
adopted by [74] then gives µ = 0.69 MPa, which is used to calculate the bifurcation pres-
sure.
It can be seen that for a very thin tube (0 < H/D < 0.1), bifurcation occurs at
a small external pressure. For tubes with larger wall thickness, when H/D > 0.1, the
bifurcation pressure increases rapidly. For 0 < H/D < 0.4, our results are in accord with
the experimental results of [74] and von Mises’ formula. When H/D > 0.4, the latter
curve increases more rapidly than our results.
Although we have considered a tube of finite length, a limitation of the present study
is that we have initiated the bifurcation analysis from a deformed circular cylindrical con-
figuration and adopted rather special incremental boundary conditions on the ends of the
tube. These might prevent realistic post-buckling behaviour for which large deformations
can occur in either the axial or azimuthal direction near the ends. Thus, our results only
apply for the initial bifurcation behaviour. Many interesting phenomena, such as self-
exited oscillations in collapsible tubes conveying fluid, occur in the post-buckling phase,
where the cross-sectional area typically takes on an elliptical or dumbbell shape. These
are excluded in the present analysis.

4.6 Conclusion

Axisymmetric and asymmetric bifurcations of circular cylinders under external pressure

combined with axial loading have been analyzed in detail using a particular model strain-
energy function appropriate for nonlinear elastic deformations of rubberlike materials.
Unlike the models used by [72] and [77], which are applicable only for thin-walled tubes,
this study presents results for a wide range of tube wall thickness on the basis of the exact
3D theory of finite elasticity. A more general description of the bifurcation behaviour of
thick-walled tubes subject to external pressure combined with axial loading, including axial
compression and extension, has been presented. Good agreement with previous studies has
been found, and extensive comparisons with results for thin-shell theory are made. Our
results show that the critical bifurcation pressure deviates from the thin shell prediction
CHAPTER 4. BIFURCATION 59

150
L/B=5
L/B=10
Mises formula

100

P
50

see Fig. (b)

0
0.2 0.4 0.6 0.8
H/D

(a)

L/B=5
14
L/B=10
Mises formula
12

10

8
P

0
0 0.05 0.1 0.15 0.2 0.25
H/D

(b)

Figure 4.13: (a) Mode m = 2 bifurcation pressures P vs. H/D for silicone rubber tubes
for λz = 1.005; the continuous curve is for L/B = 10, and dash-dot curve is for L/B = 5.
The dashed curve corresponds to von Mises’ theoretical result. (b) the enlarged area
indicated in (a). The symbols are from the Weissman and Mockros experimental results:
∇ represents bifurcation points at 50% volume collapse and 4 at 70%.

in both the very thin and thick-walled regimes. For very short and sufficiently thick
tubes, transition from lower to higher modes occurs in the range of axial compression. We
have also shown that, contrary to thin-shell theory, for sufficiently thick tubes, transition
from lower to higher modes does not occur for sufficiently short tubes. Instead, mode 2
bifurcation becomes the sole dominant mode.
Chapter 5

Nonlinear axisymmetric
deformations

In this chapter, we restrict our attention to the nonlinear axisymmetric deformations of

elastic tubes under external pressure.

5.1 Basic equations

We consider an initially stress-free thick-walled circular cylindrical tube. In this reference

configuration the geometry of the tube is described in terms of cylindrical polar coordinates
R, Θ, Z by
A ≤ R ≤ B, 0 ≤ Θ ≤ 2π, 0 ≤ Z ≤ L, (5.1)

where A and B, respectively, are the inner and outer radii and L is the length of the
tube. Let ER , EΘ , EZ denote the associated unit basis vectors. The deformed geometry
is described in terms of cylindrical polar coordinates r, θ, z with corresponding unit basis
vectors er , eθ , ez . In what follows we shall consider axisymmetric deformations of the tube.

5.1.1 Deformation

Let u denote the displacement vector, which, for axisymmetric deformations, may be
expressed in the form
u = u(R, Z)er + w(R, Z)ez . (5.2)

60
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 61

The deformation gradient tensor F is given by

u
F = (1 + uR )er ⊗ ER + uZ er ⊗ EZ + (1 + )eθ ⊗ EΘ
R
+ wR ez ⊗ ER + (1 + wZ )ez ⊗ EZ , (5.3)

where the subscripts R and Z on u and w indicate the partial derivatives ∂/∂R and ∂/∂Z,
respectively. The matrix representation of (5.3) with respect to both sets of cylindrical
polar coordinates is  
1 + uR 0 uZ
 
 
F = 0 1 + u/R 0.
 
wR 0 1 + wZ
Using (5.3), we may calculate the right Cauchy-Green deformation tensor, defined by
C = FT F, where T denotes the transpose. This yields

C = [(1 + uR )2 + wR
2
]ER ⊗ ER + (1 + u/R)2 EΘ ⊗ EΘ + [u2Z + (1 + wZ )2 ]EZ ⊗ EZ

+ [uZ (1 + uR ) + (1 + wZ )wR ](ER ⊗ EZ + EZ ⊗ ER ). (5.4)

We also note the polar decomposition (2.6) discussed in Chapter 1, where R is a proper
orthogonal tensor and U is the right stretch tensor, which is positive definite and symmet-
ric. Thus, C = U2 . The eigenvalues of U are the principal stretches of the deformation,
denoted λi , i = 1, 2, 3. The principal axes of C and U coincide and we can see immediately
from (5.4) that EΘ is a (Lagrangian) principal axis, which corresponds to the principal
stretch λ2 = 1 + u/R. The other two principal axes lie parallel to the (R, Z) plane and
can be defined in terms of an angle ψ via

E0R = cos ψER + sin ψEZ , E0Z = − sin ψER + cos ψEZ . (5.5)

The connection between principal and reference axes can be seen clearly in Fig. 5.1.
The corresponding principal stretches are taken as λ1 and λ3 , respectively. Then, we
have
C = λ21 E0R ⊗ E0R + λ22 EΘ ⊗ EΘ + λ23 E0Z ⊗ E0Z . (5.6)

5.1.2 Material properties and equilibrium

The material of the tube is considered to be incompressible, so that the constraint

J = det F = det U = λ1 λ2 λ3 ≡ 1 (5.7)

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 62

EZ
0
EZ

0
0
EΘ = EΘ ER

ψ
ER

Figure 5.1: The connection between principal and reference axes.

must be satisfied for every material point X. Subject to this constraint, the elastic proper-
ties of the material can be described in terms of a strain-energy function W (F), defined per
unit volume. By objectivity W (F) = W (U). The associated Biot stress tensor, denoted
here by T, is then given by
∂W
T= − pU−1 , (5.8)
∂U
where p is a Lagrange multiplier associated with the constraint (5.7). For details of the
Biot stress tensor we refer to [60]. For the considered deformation p is a function only of
R and Z.
Now, for an isotropic material W is a function only of the principal stretches λ1 , λ2 , λ3 ,
again subject to (5.7), and T has the same principal axes as U. The principal Biot stresses
are then simply
∂W
ti = − pλ−1
i , i = 1, 2, 3. (5.9)
∂λi
Let S denote the nominal stress tensor. Then, since the material is isotropic, we have

S = TRT , (5.10)

where R is obtained from the polar decomposition as R = FU−1 . As mentioned in

Chapter 1, in the absence of body forces the equilibrium equation is expressed in terms of
the nominal stress as
DivS = 0, (5.11)
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 63

where Div is the divergence operator with respect to X. Alternatively, in terms of the
Cauchy stress tensor, denoted σ and given by σ = J −1 FS, the equilibrium equation may
be written equivalently as
divσ = 0. (5.12)

The principal Cauchy stresses are given by

∂W
σi = λi − p, i = 1, 2, 3. (5.13)
∂λi
On the external lateral surface of the tube a pressure P , per unit deformed area, is
applied, while the inner surface is kept free of traction. The boundary conditions on these
surfaces may then be given as

 −P F−T N on R = B
ST N = (5.14)

0 on R = A,

where N is the unit outward normal to the lateral surface of the tube in the reference
configuration, i.e. N = ER on R = B and N = −ER on R = A.
On the ends of the tube the displacement is taken to vanish except for the special case
in which we consider the deformation to maintain circular symmetry. Thus,

u = w = 0 on Z = 0, L. (5.15)

For the linear and nonlinear cases, the boundary conditions are illustrated in Fig. 5.2.
For the specific calculations we make use of the neo-Hookean strain-energy function,
which is given by (2.54).

5.2 Linear and nonlinear equations

We consider the nonlinear formulation with the boundary conditions specified above to-
gether with two special cases: the first is nonlinear but assumes that the deformation is
radially symmetric, for which an analytical solution is obtained, while the second is based
on the linear theory of elasticity. These special cases serve to verify our C++ code and to
highlight, in particular, the differences between the linear and nonlinear results.

5.2.1 Radially symmetric case

If the deformation is radially symmetric then the deformed geometry has the form

a ≤ r ≤ b, 0 ≤ θ ≤ 2π, 0 ≤ z ≤ l, (5.16)
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 64

u=0; w=0

P n P
N
dA da

(1) Reference configuration (2) Current configuration

Figure 5.2: Boundary conditions for linear and nonlinear cases (a) in Reference configu-
ration (b) in current configuration.

where a and b, respectively, are the deformed inner and outer radii of the tube and l is its
length.
For this special case, we assume that the displacement is given by u = u(R)er , so that
there is no dependence on Z and w is identically zero. Then the deformation gradient
tensor F in (5.3) specializes accordingly, and the right Cauchy-Green deformation tensor
in (5.4) reduces to

C = (1 + uR )2 ER ⊗ ER + (1 + u/R)2 EΘ ⊗ EΘ + EZ ⊗ EZ . (5.17)

It follows that the Lagrangian principal axes coincide with the basis vectors ER , EΘ , EZ
and the principal stretches are
u
λ1 = 1 + uR , λ2 = 1 + , λ3 = 1. (5.18)
R
Furthermore, S = T and hence

S = t1 ER ⊗ ER + t2 EΘ ⊗ EΘ + t3 EZ ⊗ EZ . (5.19)

The equilibrium equation (5.11) specializes to the single component

1
SRr,R + (SRr − SΘθ ) = 0, (5.20)
R
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 65

where SRr = t1 , SΘθ = t2 and ,R ≡ d/dR. For the neo-Hookean material (2.54) we then
obtain, on use of (5.18),

t1 = µλ1 − pλ2 , t2 = µλ2 − pλ1 , t3 = µ − p, (5.21)

where the incompressibility condition λ1 λ2 = 1, or equivalently

u + (R + u)uR = 0, (5.22)

has been used. The latter can be integrated to give r = R + u in the form

r2 = R2 + a2 − A2 . (5.23)

The component form of the boundary condition (5.14) may now be written

 −P λ2 on R = B
SRr ≡ t1 = (5.24)

0 on R = A.
Using (5.21) and noting that Rλ2,R = λ1 − λ2 we may integrate (5.20) and use the
boundary conditions (5.24) to obtain
µ ¶ µ 2 ¶
Ab 1 A B2
P = µ ln + µ − 2 . (5.25)
Ba 2 a2 b

5.2.2 The linear case

In the linear theory of incompressible isotropic elasticity the (Cauchy) stress tensor is given
by

σ = −pI + µ[gradu + (gradu)T ], (5.26)

where I is the identity tensor.

Then, for the axisymmetric situation, the equilibrium equation (5.12) has two compo-
nents that are not satisfied trivially, namely
1
σrr,r + σzr,z + (σrr − σθθ ) = 0, (5.27)
r
1
σrz,r + σzz,z + σrz = 0, (5.28)
r
and the incompressibility constraint is
u
ur + + wz = 0. (5.29)
r
The boundary conditions (5.15) are unchanged, but (5.14) may be reduced to

 −P n on R = B
σn = (5.30)

0 on R = A,
there being no distinction between the deformed and reference configurations.
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 66

5.2.3 The nonlinear case

Comparing (5.4) and (5.6) and using (5.5), we obtain

λ21 cos ψ 2 + λ23 sin ψ 2 = (1 + uR )2 + wR

2
, (5.31)

λ21 sin ψ 2 + λ23 cos ψ 2 = u2Z + (1 + wZ )2 , (5.32)

(λ21 − λ23 ) sin ψ cos ψ = uZ (1 + uR ) + wR (1 + wZ ), (5.33)

and λ2 = 1 + u/R. From (5.31)–(5.33), it follows that

(λ21 − λ23 ) cos 2ψ = wR

2
− u2Z + (1 + uR )2 − (1 + wZ )2 , (5.34)

(λ21 − λ23 ) sin 2ψ = 2[uZ (1 + uR ) + wR (1 + wZ )]. (5.35)

It turns out that we must take λ1 > λ3 , and hence we obtain

p p
2λ1 = (uZ − wR )2 + (uR + wZ + 2)2 + (uZ + wR )2 + (uR − wZ )2 , (5.36)
p p
2λ3 = (uZ − wR )2 + (uR + wZ + 2)2 − (uZ + wR )2 + (uR − wZ )2 . (5.37)

Recalling that the Biot stress tensor has the same principal axes as U we may write

T = t1 E0R ⊗ E0R + t2 EΘ ⊗ EΘ + t3 E0Z ⊗ E0Z , (5.38)

and hence from (5.10) with R = FU−1 , we obtain the components of the nominal stress
tensor in the form

SRr = (λ−1 −1 −1 2 −1 2
1 t1 − λ3 t3 )uZ sin ψ cos ψ + (1 + uR )(λ1 t1 cos ψ + λ3 t3 sin ψ),

SRz = (λ−1 −1 −1 2 −1 2
1 t1 − λ3 t3 )(1 + wZ ) sin ψ cos ψ + wR (λ1 t1 cos ψ + λ3 t3 sin ψ),

SZr = (λ−1 −1 −1 2 −1 2
1 t1 − λ3 t3 )(1 + uR ) sin ψ cos ψ + uZ (λ1 t1 sin ψ + λ3 t3 cos ψ),

SZz = (λ−1 −1 −1 2 −1 2
1 t1 − λ3 t3 )wR sin ψ cos ψ + (1 + wZ )(λ1 t1 sin ψ + λ3 t3 cos ψ), (5.39)

together with SΘθ = t2 .

The appropriate specialization of the equilibrium equation (5.11) then yields the two
equations
1
SRr,R + SZr,Z + (SRr − SΘθ ) = 0, (5.40)
R
1
SRz,R + SZz,Z + SRz = 0, (5.41)
R
the incompressible condition is

(1 + u/R)[(1 + uR )(1 + wZ ) − uZ wR ] = 1, (5.42)

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 67

and the boundary condition (5.14) specializes to


 −P (1 + u/R)(1 + wZ ) on R = B
SRr = (5.43)

0 on R = A,

with

 P (1 + u/R)uZ on R = B
SRz = (5.44)

0 on R = A.

5.3 Finite element algorithm

To solve the nonlinear partial differential equations, the object-oriented package Libmesh
[44] is used, which is a framework for solving and analyzing systems of nonlinear equations
using the finite element method. It is also an interface to the high quality software PETSc,
which is used to solve linear systems on both serial and parallel platforms.

5.3.1 Discretization

We discretize the governing PDEs (5.11) with the constraint (5.7) using the weighted
residual-Galerkin method. The elastic domain is divided into a set of sub-domains.
Libmesh offers the options of quadratic elements of 9-node quadrilateral and 6-node trian-
gle type. Using a mixed interpolation approach, the displacement components u, w and the
radial coordinate R are interpolated by quadratic shape functions Ni , while the Lagrange
multiplier p is interpolated by linear shape functions Li , i.e.
n1
X n1
X
u= Nk (ξ, η)uk , w= Nk (ξ, η)wk ,
k=1 k=1
n1
X n2
X
R= Nk (ξ, η)Rk , p= Lk (ξ, η)pk ,
k=1 k=1

where n1 , n2 are the element node numbers, which are dependent on the element type
chosen, and ξ and η are natural coordinate variables, corresponding to isoparametric finite
elements.
This allows us to write the discretized nonlinear governing equations as

< = K(U)U − F(U) = 0, (5.45)

where U is the global vector of unknowns, K(U) is the global stiffness matrix, F(U)
denotes the force vector, which is also dependent on U, and < is the global residual
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 68

vector, which should be 0 for an exact solution. Note that U was used earlier for the right
stretch tensor, which does not appear hereon so there is no conflict of notation. Numerical
simulations show that the 6-node triangle is more efficient than the 9-node quadrilateral
element for large distortions. The formulation of the finite element matrices is problem
dependent, as shown in Section 4.3 below.

5.3.2 Newton’s method

To solve systems of nonlinear equations such as (5.45), the SNES library of PETSc [4] is
called by Libmesh. The SNES library provides a powerful suite of numerical routines, and
Newton-Krylov methods provide the core of the package, including line search and trust
region techniques. Newton’s iteration may be implemented by

Ur+1 = Ur − J−1 (Ur )<(Ur ), (5.46)

where r is the iteration number and J is the Jacobian matrix, which, by using (5.45), is
defined by
∂<(Ur ) ∂K(Ur ) ∂F(Ur )
J(Ur ) = = K(Ur ) + Ur − . (5.47)
∂U ∂U ∂U
Convergence is achieved when the relative residual tolerance ||<(Ur )||/||<(U0 )|| (in the
l2 norm) is less than 10−8 or the absolute tolerance ||<(Ur )|| is less than 10−12 , where
<(U0 ) is the initial residual.

5.3.3 Detailed discretizing integrations

Radially symmetric case

Applying Galerkin’s method to equation (5.20), we obtain

Z Z
1
Ni SRr,R dΩ + Ni (SRr − SΘθ )dΩ = 0, (5.48)
Ω Ω R
where Ω is the integration domain. The domain of integration Ω is the physical domain
in the reference configuration corresponding to the (R, Z) tube section. For each element,
(5.48) can be integrated by parts to give
Xn1 Z Z
1
− ( Ni Nj + RNi,R Nj,R )dRdZuj
j=1R Z R
n2 Z Z
X n1
X n1
X
+ Lj [(1 + Nk,R uk )Ni + (R + Nk uk )Ni,R ]dRdZpj
j=1 R Z k=1 k=1
Z Z Z
=− (RNi SRr )|R
R1 dZ
2
+ (Ni + RNi,R )dRdZ. (5.49)
Z R Z
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 69

Equation (5.22) may be discretized similarly to give

n1 Z Z
X n1
X
RLi [Nj + (R + Nk uk )Nj,R ]dRdZuj = 0. (5.50)
j=1 R Z k=1

Here we have adopted 9-node quadrilateral elements in order to achieve better accuracy,
so that n1 = 9, n2 = 4. Note that when assembled globally the boundary integrals in
(5.49) cancel out except at the boundaries of the tube.

The linear case

The discretized equations for the linear case can be obtained by using a similar procedure
to that used for the radially symmetric case. This yields
n1 Z Z
X X1 n Z Z
µ
µr(2Ni,r Nj,r + Ni,z Nj,z ) + 2 Ni Nj drdzuj + µrNi,z Nj,r drdzwj
r z r r z
j=1 j=1
n2 Z
X Z Z Z
− (rNi,r + Ni )Lj drdzpj = (rNi σrr )|rr21 dz + r(Ni σzr )|zz21 dr,(5.51)
j=1 r z z r

n1 Z Z
X n1 Z Z
X
µrNi,r Nj,z drdzuj + µr(Ni,r Nj,r + 2Ni,z Nj,z )drdzwj
j=1 r z j=1 r z
n2 Z Z
X Z Z
− rNi,z Lj drdzpj = (rNi σrz )|rr21 dz + r(Ni σzz )|zz21 dr, (5.52)
j=1 r z z r

n1 Z Z
X X1 n Z Z
1
rLi (Nj,r + Nj )drdzuj + rLi Nj,z drdzwj = 0. (5.53)
r z r r z
j=1 j=1

For the linear case, the 6-node triangle is used, so that n1 = 6, n2 = 3. This is also used
for the following nonlinear case since for large distortions the triangular element shows its
superiority over the rectangular element.

The nonlinear case

On applying Galerkin’s method to equations (5.40)–(5.42) and integrating by parts, we ob-

tain the stiffness matrix, which can be written in many different ways since the dependent
variables are nonlinearly coupled in each of the terms of the stiffness matrix. In general,
in the nonlinear case we obtain the discretization by separating off the terms that also
appear in the radially-symmetric and linear cases so that the final equations in the non-
linear case can be taken as the corresponding linear ones multiplied by some complicated
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 70

higher order coefficients. The final forms of the discretized equilibrium equations and the
incompressibility condition are

Xn1 Z Z ½· ¸ ¾
µ (1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
Ni Nj + RNi,R + Nj,R + sin 2ψ( − )Nj,Z
R 2λ1 2λ3 2 λ1 λ3
j=1 R Z
½ · ¸ ¾
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
+ RNi,Z sin 2ψ( − )Nj,R + + Nj,Z dRdZuj
2 λ1 λ3 2λ1 2λ3
Xn2 Z Z · ¸ Z Z
1 1 + cos 2ψ −2 (5.54)
− R Ni + λ1 Ni,R Lj dRdZpj = − µNi dRdZ
R Z R+u 2 R Z
j=1
Z Z · ¸ Z Z
1 + cos 2ψ (1 − cos 2ψ)t3 1 t1 t3
− RNi,R µ + dRdZ − RNi,Z sin 2ψ( − )dRdZ
ZR Z 2Z 2λ3 R Z 2 λ1 λ3
R2 Z2
+ (RNi SRr )|R1 dZ + R(Ni SZr )|Z1 dR,
Z R

n1 Z Z
X ½· ¸ ¾
(1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
RNi,R + Nj,R + sin 2ψ( − )Nj,Z
2λ1 2λ3 2 λ1 λ3
j=1 R Z
½ · ¸ ¾
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
+ RNi,Z sin 2ψ( − )Nj,R + + Nj,Z dRdZwj
2 λ1 λ3 2λ1 2λ3
n2 Z Z Z Z (5.55)
X 1 + cos 2ψ −2 1 t1 t3
− R λ3 Ni,Z Lj dRdZpj = − RNi,R sin 2ψ( − )dRdZ
2 2 λ1 λ3
j=1 R Z R Z
Z Z · ¸ Z Z
(1 − cos 2ψ)t1 1 + cos 2ψ R2
− RNi,Z +µ dRdZ + (RNi SRz )|R1 dZ + R(Ni SZz )|Z
Z1 dR,
2

R Z 2λ1 2 Z R

n1 Z Z
X µ ¶
1
RLi Nj + Nj,R dRdZuj
R Z R+u
j=1
n1 Z Z
X
+ RLi [(1 + uR )Nj,Z − uZ Nj,R ] dRdZwj = 0. (5.56)
j=1 R Z

Using equations (5.54)–(5.56), we obtain the stiffness matrix K. It is worth mentioning

that in order to achieve convergence of the solutions J needs to be computed analytically
from (5.47). Although a much simpler approach to estimating the true Jacobian matrix J is
commonly used by using the stiffness matrix K this does not work for our nonlinear model.
This indicates that the second and third terms in the expression (5.47) are important and
cannot be neglected.

5.4 Numerical algorithms

The work for the numerical processes are summarized as follows:

1. write the file mesh.xda for symmetric mesh generation.

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 71

2. write the main programme to call the linear/nonlinear solvers to solve the equation
systems by applying the external pressure as a sequence of increments.

3. write the subfunction void compute_jacobian to evaluate the true global Jacobian
matrix J.

4. write the subfunction void compute_residual to evaluate the global residual vector
<.

5. write the subfunction struct stress_vector cauchy_stress to evaluate the stresses

and principal stresses.

6. write the subfunction std::vector<Real>& stretch to evaluate the principal stretches.

The algorithms for nonlinear case including the main programme and a sub-programme
for constructing the global Jacobian matrix J are presented by the flowcharts as follows.
The algorithms for the main programme are illustrated in Figs. 5.3–5.5 and the algorithms
for constructing the global Jacobian matrix are shown in Figs. 5.6–5.10. We don’t provide
the algorithms for the radially symmetric and linear case since they are similar and simpler
than those for the nonlinear case.

5.5 Numerical results

To demonstrate the differences between the nonlinear and linear cases, three options will be
considered for the tube geometry: thick-walled short tubes with A/B = 0.5 and L/B = 1,
thick-walled longer tubes with A/B = 0.5 and L/B = 5, and thin-walled longer tubes
with A/B = 0.8 and L/B = 5.
Henceforth, all the variables are used in dimensionless form, but without change of
notation. The radial coordinates R and r and the displacement components u and w are
non-dimensionalized with B; the axial coordinates Z and z with L; the pressure P and
the stress components σij with the shear modulus µ.

5.5.1 Thick-walled short tubes: A/B = 0.5 and L/B = 1

Displacements and stretches

As both the linear and nonlinear models should agree when deformation is small, to
validate the numerical approach, a comparison of the nonlinear and linear models for
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 72

Main Program

Initialize libmesh

Read in control parameters from

command lines and file nonlinear.in

Read in mesh (TRI3) from mesh

file “mesh.xda”

Convert a mesh with linear elements (TRI3) into

a mesh with second-order elements (TRI6)

Create a nonlinear implicit system namely “nonlinear”

and add the variables (u, w, p) to the system

Give the system a pointer to the sub-functions

that update the residual R and Jacobian K

Give the system a pointer to the

initialization function

Initialize equation system

data structure

Figure 5.3: Flowchart for for nonlinear case, part 1.

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 73

Reset flags of both linear and nonlinear solvers

Initialize the data structures for the

equation system

Preparations for post processing

Nonlinear pressure loops

Refine the coarse mesh and

reinitialize the equation system

Assemble and solve the "nonlinear" system

Is the
post_processing
No
true?

Yes

Figure 5.4: Flowchart for for nonlinear case, part 2.

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 74

Call the sub-function:

cauchy_stress(equation_systems,*nonlin
ear_system.solution, p_step)
to get the stresses and principle stresses

Call the sub-function:

stretch(equation_systems,*nonlinear_sys
tem.solution, p_step)
to get the principle stretches

get the displacements, stresses, principle

stresses and principle stretches at a
particular point at each pressure step

Increase pressure: pressure

Output the displacements into GMV or

Tecplot format

End the pressure loops

End Main Program

Figure 5.5: Flowchart for for nonlinear case, part 3.

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 75

Matrix assembly function

void compute_jacobian (const
NumericVector<Number>& soln,

SparseMatrix<Number>& jacobian)

Get a reference to the equation system

Get pressure

Get const reference to the

mesh

Get a reference to the

"nonlinear" system

Get numerical ids for each

variable

Get the Finite Element

types for variables

Build finite element objects of the

unknown variables and a special finite
element object for boundary
integration

Figure 5.6: Flowchart for the Global stiffness matrix, part 1.

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 76

Define the Gauss quadrature rule

Define Jacobian*quadrature weight :JxW,

element shape functions and shape
function gradients.

Define the data structures to contain

element matrices

Define vectors to contain the DOF

for the element and variables

Loop over all the elements in the mesh.

(This involves compute element matrix and right-
hand-side contribution.)

Fill the DOF vectors with global degree of

freedom indices for the element.

Figure 5.7: Flowchart for the Global stiffness matrix, part 2.

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 77

Compute element data for the current

element. This involves computing the location
of quadrature points and the shape functions
(phi, dphi)

Zero the element matrices before

summing them.

Reposition the sub-matrices

Quadrature point loops

Define values to hold the solution & its

gradient at the previous pressure iterate and
Compute the variables and their gradient from
the previous pressure iterate.

Compute sub-matrices for

current element

Figure 5.8: Flowchart for the Global stiffness matrix, part 3.

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 78

End quadrature point loops

Impose boundary condition

Define penalty value

Dose the current No

element has no Nothing
neighbour on a side? need do

Yes

Impose top and bottom boundary

condition using penalty method

Define shape function (phi_face,

dphi_face) and
Jacobin*quadrature weight
(JxW_face) on the face of the
element

Boundary id==1?
(1=right side of element)

Yes

Figure 5.9: Flowchart for the Global stiffness matrix, part 4.

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 79

Quadrature point loops on

the face of the element

Define and compute values to hold solution

and its gradient at previous pressure iteration

Modify the element matrix

contribution

End quadrature point loops on the

face of the element

End boundary condition

Add the element matrix to the global matrix

End of element loop

End of program

Figure 5.10: Flowchart for the Global stiffness matrix, part 5.

CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 80

small pressure (P = 0.05) has been made, as shown in Fig. 5.11, in which contour plots
of the values of u and w for each case are illustrated, superimposed on the deformed (r, z)
section 1 of the tube. As expected, the distributions of the displacement components u
and w for these two cases are virtually indistinguishable. However, the difference between
the nonlinear and linear models increases as the pressure increases. This is highlighted in
Fig. 5.12(a), where the displacement u in the radial direction versus the external loading
P at point (R, Z) = (0.5, 0.5) is shown.
Figure 5.12(a) shows that the linear theory overestimates the displacement u, espe-
cially for large external pressure (P & 1.5). For example, for P = 2, the predictions of
u for the linear and nonlinear cases are 0.416 and 0.291, respectively, an overestimate of
43%. Further validation of our numerical code can be made by comparing the analyti-
cal and numerical solutions for the radially-symmetric case shown in Fig. 5.12(a). The
curves are indistinguishable in this figure. Note that the radially-symmetric and nonlin-
ear curves intersect at P ≈ 1.15. For P & 1.15, u increases with P more slowly for the
radially-symmetric case than for the nonlinear one, thus significantly underestimating the
displacement in the radial direction. The linear theory predicts a smaller axial displace-
ment w than the nonlinear theory, while for the radially-symmetric case w = 0; see Fig.
5.12(b). The differences in the results for the considered point are representative of those
seen at other points, details for which are not shown here.
The deformation, as distinct from the displacement, can be characterized in terms of
the stretches, and this is illustrated in Fig. 5.13, which shows how the principal stretches at
the point R = 0.5, Z = 0.5 vary with the pressure P . It can be seen that at this point λ1 >
λ3 > 1 and λ2 < 1. Again, for smaller pressure (P . 0.4), the principal stretches are almost
the same for the linear and nonlinear cases. It is clear, and of course not surprising, that the
linear theory provides an accurate prediction only for small deformations, corresponding
to the maximum principal stretch λ1 less than about 1.1. However, as we shall see in
the next section, the linear–nonlinear correspondence reduces to λ1 less that about 1.05,
i.e. to a strain of about 5%, when the stress components are considered. As the pressure
increases the nonlinear theory predicts larger values of λ1 , λ2 and λ3 than the linear theory.
It should be remarked that the incompressibility condition λ1 λ2 λ3 = 1 is violated for the
principal stretches calculated for the linear theory except for very small values of P . This
just emphasizes that the linear theory cannot be expected to be valid except for small
1
Note that the displacements u and w are too small to be seen in Fig. 5.11.
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 81

1 1

0.9 0.9

0.8 0.8
-0.000651856 -0.000644825
-0.00130371 -0.00128965
0.7 -0.00195557 0.7 -0.00193447
-0.00260742 -0.0025793
-0.00325928 -0.00322412
0.6 -0.00391114
0.6
-0.00386895
z 0.5
-0.00456299
-0.00521485
z 0.5
-0.00451377
-0.0051586
-0.00586671 -0.00580342
-0.00651856 -0.00644825
0.4 -0.00717042 0.4
-0.00709307
-0.00782227 -0.0077379
-0.00847413 0.3 -0.00838272
0.3 -0.00912599 -0.00902755
-0.00977784 -0.00967237
0.2 0.2

0.1 0.1

0 0
0.5 0.6 0.7 0.8 0.9 1 0.5 0.6 0.7 0.8 0.9 1

r r
(a) (b)

1 1

0.9 0.9

0.8 0.8
0.00186801 0.00189531
0.00160115 0.00162455
0.7 0.00133429 0.7 0.00135379
0.00106744 0.00108304
0.000800576 0.000812276
0.6 0.000533718
0.6
0.000541518
z 0.5
0.000266859
0
z 0.5
0.000270759
0
-0.000266859 -0.000270759
-0.000533718 -0.000541518
0.4 -0.000800576 0.4
-0.000812276
-0.00106744 -0.00108304
-0.00133429 0.3 -0.00135379
0.3 -0.00160115 -0.00162455
-0.00186801 -0.00189531
0.2 0.2

0.1 0.1

0 0
0.5 0.6 0.7 0.8 0.9 1 0.5 0.6 0.7 0.8 0.9 1

r r

Figure 5.11: Deformed profiles of the tube section in (r, z) space for a tube with A/B =
0.5, L/B = 1 subject to external pressure P = 0.05 with the displacement distributions
u and w superimposed: (a) linear u; (b) nonlinear u; (c) linear w; (d) nonlinear w. The
plots correspond to R ∈ [0.5, 1], Z ∈ [0, 1].

pressures and the accompanying small deformations.

To better understand the effect of the nonlinear contributions in equations (5.40)–
(5.42), the displacement distributions are plotted for a relatively large value 2.3 of the
pressure P in Fig. 5.14. For the purpose of comparison, the corresponding linear results
are also shown. Some significant differences between the linear and nonlinear models can
be observed in Fig. 5.14. The displacement in the radial direction is so large in the linear
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 82

0
0.16

0.14

-0.1
0.12 nonlinear
0.1

u -0.2 w 0.08
radially symmetric
0.06 linear
-0.3
0.04

nonlinear 0.02
-0.4
linear
0
0 1 2 3 4 5 0 0.5 1 1.5 2 2.5
P P

(a) (b)

Figure 5.12: Plots of displacement against pressure for a tube with A/B = 0.5, L/B = 1 at
specific points: (a) u versus P at point (R, Z) = (0.5, 0.5); radially symmetric (dash-dotted
curve); linear (dashed line); nonlinear (solid curve): (b) w versus P at point (R, Z) =
(0.5, 0.75); linear (dashed line); nonlinear (solid curve).

1.8
1

1.6
3
1

1.4 3

1.2

λ 1

0.8

0.6
2
0.4 2

0.2
0 0.5 1 1.5 2 2.5
P

Figure 5.13: Plots of λ1 , λ2 and λ3 (labelled 1, 2, 3, respectively) versus P calculated at

the point (R, Z) = (0.5, 0.5) for a tube with A/B = 0.5, L/B = 1: linear (dashed lines);
nonlinear (solid curves).
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 83

case that the middle section of the tube almost comes into self contact on the axis R = 0.
For the nonlinear case, the most striking feature is the bulging out at the corners, which
is barely visible in the linear case. This causes the displacement pattern and magnitude
to change. The radial displacement u changes between −0.47 and 0 in the linear case,
and between −0.31 and 0.0145 in nonlinear case. The axial displacement w has the range
−0.095 to 0.095 (linear case) and −0.15 to 0.15 (nonlinear case). This is consistent with
the corner bulging at R = 0.5 on the ends of the tube, which stretches the tube section in
two opposite axial directions.
-0
-0 6 3

-0. -0.0
1 1
.19
.16

06 37 -0.0062 -0 -0.0062
95 8 .16 - -0.0378
46

-0.0378 46 0.13 -0.1

-0.06 29 012 -0.069
-0 -0.10 5
.2 -0 .
-0.1329 12 95
9 13 25
97 -0.2
-0 280 -0.1 -0.00615694
0.8 .3 -0.1 646 0.8
-0

2 963 -0.00615694 -0 -0.0378445

.3

30 .19
8

-0.0378445 63
64

-0.0695321

-0
-0 -0.1

.2
-0.069532 .2
-0.4

-0 64 -0.10122

5
.3 2 6

97
5 80
-0.10122 -0.132907
49

47
8

-0.29
-0.132907
-0

0.6 0.6 -0.164595

-0

.2
.2

-0.164595
5
9

-0.196282
97

13
13

-0.196282
z -0.22797 z -0.22797
-0.259658
97
30

-0.259657 -0.291345
25
.32

63
81

0.4 0.4 9
-0 .

-0.291345
-0

.1 -0.323033
64
41

80 -0
8
-0 .

22 -0.323032
.3

0.
-0.35472
-0

- 80

97
47 3 -0.35472 .2
2 -0.386408

.25
.3
5 3 96 -0
-0 91 -0.
1 -0.386407 646 -0.418095

-0
-0.2 -0.1
0.2 -0.418095 0.2 -0.449783
7
2 59 -0.449783 29
.
-0 8 0 .1646 -0 .13
-0 9
30

22 -0.132 63
-0 . 6 3 -0.0695
2

9
.3

9 -0.1012 .1 12 -0.0378
-0 10
-0

.1 -0 .
-0 9 5
-0.06 -0.0378

8
-0.0062
0 0

-0.037
-0.0062 -0.0062

0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1

r r
(a) (b)

80 789
13 0.0
1 1 0.

00
0.0000
94

0.00
0.0394
03
0.

9
0.0

8
07
97

0. 3
19

0 18
0.01

00 0.1
0.0
0.8 0.138016 0.8 0.138016
1 0.118299 0.1183
59 9 59
1
0.0985831
0.0 0.0985829 0.078 0.0
94 0.0788663 0.0788665
0.03
0.0591499
-0

0.6 0.0591497 0.6 0.0394

.01

0.0197 0.0197 0.0394332

0.0000

0.0394331
97

0.0197166 0.0197166
z 0.0000 0.0000
z 0.0000
0.0000

0 0
-0.0197166 -0.0197 -0.0197166
0.0000

0.4 -0.0 -0.0394331 0.4 -0.0 -0.0394332

19 394
7 -0.0591497 -0.0591499
-0.0
-0.0 -0.0788663 59
1
-0.0788665
59 -0.0985831
1 -0 -0.0985829
.0
3 -0.1183
0.2 94 -0.118299 0.2 -0 . -0.138016
-0.0

0.0 -0.138016 11
-0 19 83
97
78

.0 0.0 7
1

7 00
.0
9

89 0
-0
-0.
09

-0.0

-0 .
0 01 0
86

97
0.0000
39
4

0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1

r r
(c) (d)

Figure 5.14: Deformed profiles of the tube section in (r, z) space for a tube with A/B =
0.5, L/B = 1 subject to external pressure P = 2.3, with distributions of the displacement
components superimposed: (a) linear u; (b) nonlinear u; (c) linear w; (d) nonlinear w.
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 84

Cauchy stresses

The stress distributions for the linear and nonlinear cases are again almost the same for
very small external pressure, but as expected they depart significantly for large pressures.
Figures 5.15 and 5.16 show the distributions of components of Cauchy stress for P = 2.3
for the linear and nonlinear cases, respectively. Negative values of the stresses are shown
with dashed curves. In all cases, the normal stress distributions in the upper half section
have mirror symmetry with the lower half, while the distribution of shear stress is anti-
symmetric.
-2

-0.2
0.17
.72
1 -1.2 1

6
-2.169
33 1

-0.93
9 962
.16
31
-1.5

-2
-1
.70
0.9 0.9
65

1
-0.7

-2.721 -2.721
6

0.8 0.8
9
.2
-0

1.57638
1.44271
1.10818
0.7 0.7 0.847909
0.639971
69

0.253108
-2.1

0.171765
33

-0.341693
-1.2

0.6 -0.296441 0.6 -0.936493

6
.31
-0.764648
z z
-1.701

-1.53129

-3
-0.765

-1.23285
0.5 -1.70106 0.5 -2.12609
-2.7209
-2.16927
-0.296

-3.3157
0.4 -2.63747 0.4 -3.9105
-3.10568
-4.5053
-3.57389 -3.316
0.3 -4.04209
0.3 -5.1001
-2 .

-5.6949
16

-4.5103
-1 .

-6.2897
0.2 0.2
23

-4.9785 -2.721 -2.721 -6.8845

3

-1.701
-0.7

0.1 0.1
65

-0.2
96
-2.637 -2.16
9
0 0 16
3.3
0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
r r
(a) (b)

-4 .
6-96.3
2.10

70 1
1 1 33 1 63

0. .70
1 2.570

1
-3.0 70 1.

0
48 0. 8
16
1.
0.9 0.9
4
23
0.
0.8 0.8 3.50453
1
0.70
0.234

1.87876 3.03726
7
22

0.7 1.05757 0.7 2.56999

-1.406

-2 .

0.236368 2.10272
-0.58483 1.63545
0.6 -1.40603 0.6 0.234
1.16818
-0.585

z 0.5
-2.22722
-3.04842
z 0.5 0.000
0.700907
0.233636
-3.86962 -0.233636
-3

-4.69082 -0.700907
.04

0.4 -5.51201
0.4 -0.2 -1.16818
8

3 4 34
-6.33321 .2 -1.63545
-2 .

-0
22

0.3 -7.15441 0.3 -2.10272

7
-1.406

-7.97561 -0.701 -2.56999

-8.7968 -3.03726
0.2 -9.618
0.2
-3.50453
7
22
-2 .

0.1 0.1 -1.1

-0

68
-1 .

.2

-0.
3
-2.227

7
40

-3.8 -0 -2-2 01
-0.2 .701
1.

.1.5
6

70 037
0 0
87

34 0
9

0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1 0 0.1 0.2 0.3 0.4 0.5 0.6 0.7 0.8 0.9 1
r r

Figure 5.15: Deformed profiles of the tube section in (r, z) space for a tube with A/B =
0.5, L/B = 1 subject to external pressure P = 2.3, with distributions of the Cauchy stress
components superimposed. Linear case with positive contours (solid), negative contours
(dashed): (a) σ11 , (b) σ22 , (c) σ33 , (d) σ13 .
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 85

0.
-0.44832
6 -1.7725 89
1 0.40 1 19 4 -1.9 665

-3 .
-0.020 -0.8 .4
0 .
-6

15
-1.2 -1
99

-0.8

8
0.9 0.9

73

5
72
20
-2 .

-1 .
57

-0.0
3
0.8 1.68466 0.8 0.934056
1.25841

73
0.349543
0.832149

-2.5
0.7 0.7 -0.234971

7
0.405892

-0.44
5 2
.1 -0.819484

9
73
-2 -0.0203644 8

9
15

.2
-3 . -1.404

-0.8
0.6 0.6

-1
-0.446621
-1.98851
z -0.872878 z

-1.725
-2.57302
0.5 -1.29913 0.5 -3.15754
-1.72539
-3.74205
0.4 -2.15165 0.4 -4.32657
-2.57791
-4.91108

-2.5
-3.00416
0.3 0.3 -5.49559

73
-1 -3.43042

-2
.29 -6.08011

.1
9 -1.989
-0.4

-0.8 -3.85668

5
-6.66462

2
73
0.2 0.2
47

-4.28293
-7.24913

-1 .
72
5
0.1 0.1 -1.4
3 04

-7 . 5 8
87 9

7
-0 .
0. 1.29 -0.02 0
-1.7

9
.1
44
0 02 -3.00 0

24
0 - - 4 -3.74

-3
-0 .
52
25

6 -0.23
0.40
0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1
r r
(a) (b)
91

-8.507 .889

-2.3
1 1

-0.759
-3.5

9-32
89

-4.2

18
-2.8
-2.1

0.9 0.9
86

-0.308
-2 .

03
18

0.1
0.8 0.8
6

1.32526
2.56423
-0.782
-0.079

-1.48

0.622944
2.15395
0.7 -0.0793687 0.7

0
0.00
1.74368
4

-0.781682
1.3334

0
-1.48399 3
0.6 0.6

0.00
10 0.923124

00
-2.18631 0.
z z

0.0
0.512846
-2.88862
-2.889
-2.186

0.5 0.5 0.000 0.102569

-3.59093
-0.307708
-4.29325

0.0
-0.717985
0.4 0.4

00
-4.99556
-1.12826

0.00
-5.69787 0.00
-1.53854

0
-0.079

-0.782

0.3 -6.40018 0.3

-1.484

-1.94882
0

-7.1025
-2.35909
-7.80481
0.2 0.2 -2.76937
-8.50712

8
0.1 0.1

.30
84
-2.186
89

-1.4

-0
0.
-2.8

10
-4.996

-5 .
-7.805

69
3

0 -6.400 8 0
0.
92
3

0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1

r r
(c) (d)

Figure 5.16: Deformed profiles of the tube section in (r, z) space for a tube with A/B =
0.5, L/B = 1 subject to external pressure P = 2.3, with distributions of the Cauchy
stress components superimposed. Nonlinear case with positive contours (solid), negative
contours (dashed): (a) σ11 , (b) σ22 , (c) σ33 , (d) σ13 .

It can be seen from Fig. 5.15 and Fig. 5.16 that in both the linear and nonlinear cases,
the areas of stress concentration are located at the four corners. However, for the linear
case, the normal stresses σ11 , σ22 , σ33 are mostly negative, with the peak negative stresses
occurring at the two inner corners. The peak positive stresses are at the two outer corners,
and in the radial direction (i.e. for σ11 ). On the inner surface the stress σ11 is positive only
near the ends. This shows that most of the section is under compression when subject
to an external pressure. For the nonlinear case, on the other hand, all the peak stresses
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 86

(positive and negative) occur at the inner corners. This is directly due to the fact that
the inner corners are significantly squeezed out, which causes both significant compression
and tension there. The shear stress distribution σ13 is most interesting; for the linear case,
σ13 is entirely positive in the upper half and entirely negative in the lower half, with the
zero stress line at Z = 0.5. However, in the nonlinear case, because of the corner bulging,
each half section is divided into four zones between which the shear stress changes sign.
In the upper half, the innermost section is sheared upwards, while different parts of the
outermost section are subject to either positive or negative shear stress. The opposite is
true in the lower half. The general trend for short tubes is that the magnitudes of the
stresses in the nonlinear case are smaller than the corresponding linear magnitudes, with
σ11 between −4.7 and 2.11, σ22 between −7.54 and 1.52, σ33 between −8.86 and 2.02, and
σ13 between −2.87 and 2.87. These are to be compared with the linear case: σ11 from
−5.44 to 2.04, σ22 from −7.18 to 2.03, σ33 from −10 to 2.69, and σ13 from −3.74 to 3.74.
To show how the stresses change with the external pressure at a particular location, we
plot the variation of the stress components σij with the pressure at point (R = 0.75, Z =
0.75) in Fig. 5.17. Again, the differences between the results for the linear and nonlinear
models are small if P is small enough, in this case P . 0.5 for the normal stress components
and P . 0.3 for σ13 . However, significant differences are found between the linear and
nonlinear predictions as P increases, especially in the stress components σ11 and σ13 . The
nonlinear model exhibits much smaller stress magnitudes for the same applied pressure.
It is interesting that σ13 first increases rapidly with P , but reaches a maximum around
P = 1.8 and then decreases with further increase in P , as shown in Fig. 5.17(b). This is
because as the pressure increases beyond a certain level, the corners bulge out more and
more and the second left (negative) shear zone in Fig. 5.16(d) increases in size, while the
third (positive) shear zone (where the point is located) shrinks. As a result, the positive
shear stress at this point decreases for P & 1.8.
To illustrate the response of the material locally to the external forces, plots of principal
stress versus principal stretch are shown in Fig. 5.18 for the point (R = 0.75, Z = 0.75).
We note that λ1 > 1, λ2 < 1 and λ3 < 1 at the point in question. Compared with the
linear results, the nonlinear model predicts larger magnitudes of the principal stresses for
the same principal stretches. This means the stiffness of the material at this special point
becomes larger. Note that as P increases the relationship between the principal stress σ3
and the stretch λ3 loses monotonicity. It is also noted that at the point (R = 0.75, Z =
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 87

0 0.3

-0.2

-0.4
0.2
-0.6
σ11 , σ22 , σ33

σ13
-0.8
1

-1 1 0.1

-1.2 2
2
3
3
-1.4
0
0 0.5 1 1.5 2 2.5 0 0.5 1 1.5 2 2.5
P P

(a) (b)

Figure 5.17: (a) Plots of σ11 , σ22 , σ33 (labelled as 1, 2, 3, respectively) versus P for a tube
with A/B = 0.5, L/B = 1 at point (R, Z) = (0.75, 0.75). (b) Plots of σ13 versus P . Linear
(dashed lines); nonlinear (solid curves).

0.75), the angle ψ which defines the principal directions has the constant value 31.6◦ for
the linear case, while it varies in the range 29.8◦ < ψ < 32.6◦ for the nonlinear case.

0 0
0

-0.2
-0.5 -0.5

-0.4
-1 -1
σ1 σ2 σ3
-0.6
-1.5
-1.5
-0.8
-2
-2
-1
-2.5

-1.2 -2.5

-3 0.9 0.925 0.95 0.975 1

1 1.1 1.2 1.3 1.4 1.5 0.8 0.9 1

λ1 λ2 λ3
(a) (b) (c)

Figure 5.18: Plots of the principal stresses versus the corresponding principal stretches for
a tube with A/B = 0.5, L/B = 1 at the point (R, Z) = (0.75, 0.75): (a) σ1 vs λ1 ; (b) σ2
vs λ2 ; (c) σ3 vs λ3 . Linear results (dashed lines); nonlinear results (solid curves).

5.5.2 Thick-walled longer tubes: A/B = 0.5 and L/B = 5

Next, we consider a tube with the same thickness but five times longer. In this case we
find that the u and w versus P curves are similar to those for the shorter tube observed
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 88

above. The only difference is that for both the linear and nonlinear cases the deformations
of longer tubes tend to have two humps instead of one, as suggested in Fig. 5.19, i.e.
the longer tube favours mode-2 deformations for the range of the pressure applied, while
mode-1 is preferred for the shorter tubes. Figure 5.20 shows that the differences in the
stress–pressure plots between the linear and nonlinear cases are smaller than for shorter
tubes. However, the change in σ13 for the nonlinear case is interesting. As in Fig. 5.17(b),
it follows the linear curve for small P but the range of values of P for which σ13 is positive
is much smaller in this case, and it bends downwards sharply as soon as P exceeds about
0.1.

5 -0.53 5 0.698
8 91 0.547
0.0
0.095 0.243
-0.032

0.000
1.10844
4 0.981762
4 1.00204
0.855083 -0.030 0.850213
0.728404 0.698389
0.601725 0.546565
0.000
0.475046 0.394742
3 0.348367 3 0.242918
0.0910942
z 0.000
0.221688
0.0950092
z -0.0607295
0.000 -0.212553
-0.0316697
-0.158349 -0.364377
-0.285028 2 -0.516201
2 -0.668024
0.000 -0.411707
-0.538386 -0.819848
-0.665064 -0.971672
0.030
-0.791743 -1.1235
-0.918422 1
1
-1.0451
-1.17178 0.000
-0.061
-0.213
-0 .
28 -0-0.364
-0 .668 3
.8 -0.21
0 0.34
8 5 0 20
0 0.25 0.5 0.75 1 0 0.25 0.5 0.75 1

r r
(a) (b)

Figure 5.19: Distributions of the shear stress σ13 for a tube with A/B = 0.5, L/B = 5 at
pressure P = 1.01, superimposed on the deformed profile of tube section in (r, z) space:
positive values are indicated by solid contours and negative values by dashed contours.
(a) nonlinear (b) linear.

The corresponding principal stress–stretch plots are shown in Fig. 5.21. The features
of Fig. 5.21(a, b) are similar to those for the shorter tube. However, an interesting change
in the σ3 –λ3 plot is shown in Fig. 5.21(c), where an S-shaped curve is observed. This is
associated with the complicated pattern of change in the shear zones shown in Fig. 5.19(a).
The nonlinear tube tends to bulge at the two inner corners, which, when combined with
the mode-2 humps, creates a much larger negative shear zone in the upper half of the
cylinder. The smaller bulge at the corners also causes the shear stress to be split into
negative and positive regions within each half cylinder, and the negative regions emerge
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 89

0 0.05

-0.2
0
-0.4

-0.6 -0.05

σ11 , σ22 , σ33 σ13

-0.8
-0.1
-1
1

-0.15
-1.2
1

-1.4 3

2 -0.2
2 3
-1.6
0 0.5 1 1.5 2 0 0.5 1 1.5 2 2.5
P P

(a) (b)

Figure 5.20: (a) Plots of σ11 , σ22 , σ33 (labelled 1, 2, 3, respectively) versus P for a tube
with A/B = 0.5, L/B = 5 at point (R, Z) = (0.75, 4.5). (b) Plots of σ13 versus P . Linear
(dashed lines); nonlinear (solid curves).

0 0 0

-0.5
-0.5

-0.5
-1
σ1 σ2 σ3 -1

-1.5
-1.5
-1
-2

-2
-2.5

-1.5
-2.5
-3
1 1.1 1.2 1.3 0.8 0.9 1 0.98 0.985 0.99 0.995 1

λ1 λ2 λ3
(a) (b) (c)

Figure 5.21: Plots of the principal stresses versus the corresponding principal stretches for
a tube with A/B = 0.5, L/B = 5, at the point (R, Z) = (0.75, 4.5): nonlinear case. (a) σ1
vs λ1 , (b) σ2 vs λ2 , (c) σ3 vs λ3 .

and expand as the external pressure increases. The linear case shown in Fig. 5.19(b) fails
to predict the bulging at the corners at all for this case, as a result of which there is no
shear splitting zone towards the ends, although the shear zone adjacent to the boundary
region changes its sign, presumably due to the mode-2 deformation.
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 90

5.5.3 Thinner and longer tubes: A/B = 0.8 and L/B = 5

For longer and thinner tubes with A/B = 0.8 and L/B = 5, the most interesting feature is
the occurrence of higher modes (multiple humps in the deformation) in the nonlinear case.
Four modes from mode-1 to mode-4 are observed as the external pressure P increases from
0 to about 0.66, as shown in Fig. 5.22. Mode-1 occurs for 0 < P . 0.01, transitions to
mode-2 for 0.01 . P . 0.16, to mode-3 for 0.16 . P . 0.41 and mode-4 for 0.41 . P .
0.66. For larger P modes 5, 6 and 7 were obtained, although the solution for large P that
gives rise to the higher modes is more demanding on the mesh qualities. No higher modes
except mode-2 were found for the linear model.

(a) (b) (c) (d)

Figure 5.22: (Not to scale) Nonlinear modes of deformation for a tube with A/B =
0.8, L/B = 5: (a) mode-1, at P = 0.01, (b) mode-2, at P = 0.11, (c) mode-3, at P =
0.41, (d) mode-4, at P = 0.61. The contours shown are for the radial displacement u
superimposed on the (r, z) deformed profile.

Figure 5.23 shows the distributions of all the Cauchy stress components for the non-
linear case at P = 0.22. Again, there are two major differences when compared with the
corresponding linear case (not shown). One is that the nonlinear model presents a higher
mode (mode-4 in this case), where the corresponding linear case exhibits only mode-2.
The other is the shear splitting pattern in the nonlinear model, which expands from the
two ends towards the middle section. Although this is not shown here we note that the
boundary effect is more limited to near the two ends in the linear model, with the same
sign of σ13 ≥ 0 near the upper end, and σ13 ≤ 0 near the lower end. The patterns of σ22
and σ33 are also quite interesting, with the nonlinear effects more clearly focused on the
boundaries.
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 91

5 -1.256
-0.376 -0.486
-0.266 5 -1.138
-0.687 -0.687
-0.046 -0.156 -0.387
-0.046 -0.687

-0.156
-0.838

-0.838
4 4

-0 .
0.174005 0.0642304

04

3
0.0639712 -0.0860883

1
6

.9
-0
-0.0460622 -0.236407

-0.838
-0.156096 -0.386726

-0 .
-0.266129 -0.537044

76
3 3

3
-0.376162 -0.687363
z -0.486196 z -0.837682

-0.156
-0.596229 -0.988
-0.706263 -1.13832
-0.816296 -1.28864
-0.046

2 -0.92633 2 -1.43896

3
76
-1.03636 -1.58927

-0 .
-1.1464 -1.73959

-0
-1.25643 -1.88991

.8
3
-1.36646 -2.04023

8
1 1

1 3
.9
-0
-0.156

-0.687 -0.838
-0.046
-0.687
-0.156 -0.537
-0.687 6
0 -0.046
-0.596
-0.926
-1.146 -0.486 -0.266 -0.156
0 -1.289 -0.23
0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1

r r
(a) (b)

5 -2.862 5 0 1 0.0810.243
-0.34
8

-1.010 -0.47 0.114

9

-0.804 -0.146
.5

6 0.049
-0.01
-0

-0.598 -0.081
-0.804

2 0.567718
-0.39
-0.016
4 0.430918 4 0.016
0.502836
0.000
0.2251 0.437954
0.0192822 0.373072
0 -0.186536 0.30819
-0.53 -0.461
-0.392354 0.243308
0.000
3 -0.598171 3 0.178426
-0.803989 0.113544
z -1.00981 z 0.000
0.0486615
-1.21563 -0.0162205
-1.42144 -0.0811025
2 -1.62726 2 -0.145985
0.000
-1.83308 -0.210867
-0.53
0 -2.0389 -0.275749
-2.24471 -0.340631
-0.016
-2.45053 -0.405513
0.000
-0.461 -2.65635 -0.016 -0.470395
1 1
-2.86217 -0.535277
-0.187 -0.598 0.049 -0.600159
-3.06799
-3.2738
0.049 0.049
0.114 -0.016
-0.598 19 -0.081
-0 .146
0 -1.010
-2.039
-2.656 0.0 0 0.178
0.308
-0.341 -0.081
0.6 0.7 0.8 0.9 1 0.6 0.7 0.8 0.9 1

r r
(c) (d)

Figure 5.23: Cauchy stress distributions for the nonlinear case for a tube with A/B =
0.8, L/B = 5 at pressure P = 0.22 superimposed on the deformed tube section in (r, z)
space: positive values are indicated by solid contours and negative values by dashed con-
tours. (a) σ11 , (b) σ22 , (c) σ33 , (d) σ13 .

5.6 Discussion and conclusions

We have derived the general partial differential equations in Lagrangian form governing
the large axisymmetric deformations of a thick-walled tube composed of incompressible
isotropic elastic material, without any assumptions limiting the magnitude of the defor-
mation or material nonlinearity. Comparison has been made with the corresponding linear
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 92

model for tubes with different wall-thickness and length ratios.

For small deformation the linear and nonlinear models give very similar results. How-
ever, the predictions of the linear and nonlinear models are very different under large
external pressure, and the dominant nonlinear features are the corner bulging, and, for
longer tubes, the occurrence of higher modes of deformation. Note, however, that the
higher modes for longer and thinner tubes can be associated with geometrical nonlinearity
and are not features unique to material nonlinearity [38], [36]. Although we don’t dis-
tinguish the material and geometric nonlinearities in the present study, we have observed
that material nonlinearity is more important in the shorter and thicker tubes, for which
the strains computed are larger, while geometrical nonlinearity seems to dominate in the
longer and thinner tubes, for which the strains are much smaller. The Cauchy stresses,
especially the shear stress, exhibit the greatest differences between the predictions of the
linear and nonlinear theories. Shear splitting, with alternating signs of the shear stress in
different regions is a unique nonlinear feature. As a result, the end boundary constraints
have a much stronger influence on the deformation and stresses in the rest of the tube for
the nonlinear model. This is the first systematic nonlinear treatment of externally pres-
surized thick-walled elastic tubes, albeit using the simple neo-Hookean material, and the
results may have significant implications for certain physiological applications involving
soft vessels undergoing large deformation.
The nonlinear system of equations has been solved by using the C++ based finite
element package Libmesh. It should be noted that due to the complex nature of the non-
linear equations, it was extremely difficult to obtain converged solutions numerically using
approximate Newton solvers and it was necessary to derive the Jacobian matrix J analyt-
ically, and to use the corresponding linear solution as an initial solution in order to obtain
convergence. In addition, since the geometry of the tube in the reference configuration
and the boundary conditions and external pressure condition are all symmetric about the
mid-plane of the tube, a symmetric mesh needs to be used to achieve perfectly symmetric
solutions.
We have noted that the nonlinear effects for long, thin tubes are limited to a layer
p
of width (B − A)A (see [47]) near the boundaries; see Fig. 5.22. This agrees broadly
with the examples given by [47] on the behaviour of nonlinear shell-membrane materials.
No direct results can be found in [47] for a neo-Hookean cylindrical shell under external
pressure. However, qualitative comparison is possible with the results of [69] who studied
CHAPTER 5. NONLINEAR AXISYMMETRIC DEFORMATIONS 93

the deformation of a neo-Hookean cylindrical membrane under twist. They found that at
larger values of the prescribed twist, wrinkling occurs in the interior and the membrane
remains tense near the boundaries. Although these are obtained for different boundary
conditions, the nonlinear effects such as the presence of the boundary layer and multi-
modes are similar to these shown in Fig. 5.22.
Chapter 6

Three-dimensional large
deformations

In this chapter we formulate the differential equations governing the large deformation of
the cylindrical tube subject to pressure on its external lateral surface and the end condi-
tions are zero displacement. The material is assumed to be an incompressible isotropic
neo-Hookean one. The nonlinear set of equations were derived using both cylindrical and
Cartesian coordinates. For the purpose of comparison the corresponding linear equations
are also presented (using Cartesian coordinates).

6.1 Nonlinear case: cylindrical polar coordinates

Considering the special geometry of the tube, it is natural to formulate the differential
equations based on the cylindrical coordinates. This is convenient, especially, when we
impose the hydrostatic pressure on the external lateral surface of the tube.

6.1.1 Deformation gradient

Using the basic kinematic concepts described in Chapter 2 and we have the position vectors
X, x and displacement vector u (in cylindrical coordinates) in the form

x = rer + zez , (6.1)

X = RER + ZEZ , (6.2)
x = X + u, (6.3)

94
CHAPTER 6. THREE-DIMENSIONAL LARGE DEFORMATIONS 95

where displacement vector is given by

u = uer + veθ + wez . (6.4)

Note that r, θ, z are functions of R, Θ, Z.

r = r(R, Θ, Z), θ = Θ(R, Θ, Z), z = Z(R, Θ, Z), (6.5)

u = u(R, Θ, Z), v = v(R, Θ, Z), w = w(R, Θ, Z), (6.6)

and the relation between bases of cylindrical coordinates and bases of Cartesian coordinates
gives

er = cos θe1 + sin θe2 , eθ = − sin θe1 + cos θe2 , ez = e3 ,

ER = cos ΘE1 + sin ΘE2 , EΘ = − sin ΘE1 + cos ΘE2 , EZ = E3 . (6.7)

Here we have

ei = Ei , i = 1, 2, 3. (6.8)

The gradient in cylindrical coordinate system gives

∂x 1 ∂x ∂x
Gradx = ⊗ ER + ⊗ EΘ + ⊗ EZ . (6.9)
∂R R ∂Θ ∂Z

We specialize the deformation gradient (2.2) and combined with (6.9), then

F = Gradx
∂r 1 ∂r ∂r ∂θ r ∂θ
= er ⊗ ER + er ⊗ EΘ + er ⊗ EZ + r eθ ⊗ ER + eθ ⊗ EΘ
∂R R ∂Θ ∂Z ∂R R ∂Θ
∂θ ∂z 1 ∂z ∂z
+r eθ ⊗ EZ + ez ⊗ ER + ez ⊗ EΘ + ez ⊗ EZ . (6.10)
∂Z ∂R R ∂Θ ∂Z

Or by using u, v, w, we get the deformation gradient as follows

F = Gradx = Gradu + I, (6.11)

and Gradu gives

∂u ∂θ 1 ∂u v ∂θ ∂u ∂θ
Gradu = ( −v )er ⊗ ER + ( − )er ⊗ EΘ + ( −v )er ⊗ EZ
∂R ∂R R ∂Θ R ∂Θ ∂Z ∂Z
∂θ ∂v u ∂θ 1 ∂v ∂θ ∂v
+ (u + )eθ ⊗ ER + ( + )eθ ⊗ EΘ + (u + )eθ ⊗ EZ
∂R ∂R R ∂Θ R ∂Θ ∂Z ∂Z
∂w 1 ∂w ∂w
+ ez ⊗ ER + ez ⊗ EΘ + ez ⊗ EZ . (6.12)
∂R R ∂Θ ∂Z
CHAPTER 6. THREE-DIMENSIONAL LARGE DEFORMATIONS 96

¿From (6.22), we have connections between the coordinates in current and reference con-
figurations in the form


 r cos θ = R cos Θ + u cos θ − v sin θ,


r sin θ = R sin Θ + u sin θ + v cos θ, (6.13)



 z = Z + w.

¿From the first two equations in (6.13), we obtain the important relation

r−u
θ = Θ ± arccos + 2kπ, k = 1, 2, 3... (6.14)
R
∂θ ∂θ ∂θ
Then we could easily get ∂R , ∂Θ , ∂Z , but we don’t want to express these derivatives ex-
plicitly due to the complexity of the expressions.

6.1.2 Nominal stress and Cauchy stress

The strain-energy function for neo-Hookean material is

1
W (I1 ) = µ(I1 − 3). (6.15)
2

Using the definition (2.53) the nominal stress for incompressible elastic material is spe-
cialized accordingly

∂W
S= − pF−1 = µFT − pF−1 . (6.16)
∂F

For this case, the Cauchy stress tensor could be written

∂W
σ=F − pI = µFFT − pI. (6.17)
∂F

The inverse of the deformation gradient F is given by

 

 r ∂θ ∂z ∂θ ∂z 1 ∂r ∂z ∂r ∂z r ∂r ∂θ ∂r ∂θ  

 R ( ∂Θ ∂Z − ∂Z ∂Θ ) − R ( ∂Θ ∂Z − ∂Z ∂Θ ) R ( ∂Θ ∂Z − ∂Z ∂Θ ) 


 

F−1 = −r( ∂R ∂θ ∂z ∂θ ∂z
∂Z − ∂Z ∂R )
∂r ∂z ∂r ∂z
∂R ∂Z − ∂Z ∂R
∂r ∂θ
−r( ∂R ∂Z −
∂r ∂θ
∂Z ∂R ) 
(6.18)

 

 


 r ( ∂θ ∂z − ∂θ ∂z ) − 1 ( ∂r ∂z − ∂r ∂z ) r ( ∂r ∂θ − ∂θ 
R ∂R ∂Θ ∂Θ ∂R R ∂R ∂Θ ∂Θ ∂R R ∂R ∂Θ
∂r
∂Θ ∂R)

So, the nominal stress is given by in component form

   

 ∂r
r ∂θ ∂z 
 
 r ∂θ ∂z ∂θ ∂z 1 ∂r ∂z ∂r ∂z r ∂r ∂θ ∂r ∂θ  

 ∂R ∂R ∂R 
 
 R ( ∂Θ ∂Z − ∂Z ∂Θ ) − R ( ∂Θ ∂Z − ∂Z ∂Θ ) R ( ∂Θ ∂Z − ∂Z ∂Θ ) 


 
 
 

S = µ R1 ∂Θ ∂r r ∂θ 1 ∂z − p ∂θ ∂z
−r( ∂R ∂θ ∂z
∂Z − ∂Z ∂R )
∂r ∂z ∂r ∂z
∂R ∂Z − ∂Z ∂R
∂r ∂θ
−r( ∂R ∂Z −
∂r ∂θ
∂Z ∂R ) 
.

 R ∂Θ R ∂Θ 
 
 

 
 
 


 ∂r ∂θ ∂z 
 
 r ( ∂θ ∂z − ∂θ ∂z ) − 1 ( ∂r ∂z − ∂r ∂z ) r ( ∂r ∂θ − ∂r ∂θ  
∂Z r ∂Z ∂Z R ∂R ∂Θ ∂Θ ∂R R ∂R ∂Θ ∂Θ ∂R R ∂R ∂Θ ∂Θ ∂R )
CHAPTER 6. THREE-DIMENSIONAL LARGE DEFORMATIONS 97

6.1.3 Equilibrium equations

Using the Lagrangean equilibrium equation (2.18) we have

∂Sij ∂Ei ∂ej

DivS = Ek · Ei ⊗ ej + Sij Ek · ⊗ ej + Sij Ek · Ei ⊗
∂Xk Xk Xk
∂Sij ∂Ei ∂ej
= ej + Sij Ek · ej + Sij = 0 i, j = 1, 2, 3. (6.19)
∂Xi ∂Xk ∂Xi

Non-zero components of Ek · Ei,k in cylindrical coordinate system are

1 1
E2 · E1,2 = EΘ · ER,Θ = . (6.20)
R R

The final equilibrium equations could be written as

∂S11 1 ∂S21 ∂S31 1 ∂θ 1 ∂θ ∂θ

+ + + S11 − (S12 + S22 + S32 ) = 0,
∂R R ∂Θ ∂Z R ∂R R ∂Θ ∂Z
∂S12 1 ∂S22 ∂S32 1 ∂θ 1 ∂θ ∂θ
+ + + S12 + (S11 + S21 + S31 ) = 0,
∂R R ∂Θ ∂Z R ∂R R ∂Θ ∂Z
∂S13 1 ∂S23 ∂S33 1
+ + + S13 = 0. (6.21)
∂R R ∂Θ ∂Z R

Considering the complicated connections in (6.13), (6.14) and the components of nominal
stress tensor, we can easily imagine the final form of the equilibrium equations (6.21) would
be very lengthy, which will add much difficulty in the procedure of numerical simulations.

6.2 Nonlinear case: cartesian coordinates

In order to avoid the complexity of the equilibrium equations (6.21), we obtain the corre-
sponding equation system in Cartesian coordinates as follows.

6.2.1 Deformation gradient

The position vectors of an arbitral particle of the body are denoted by

x = xi ei , X = Xi Ei , x = X + u. (6.22)

where i = 1, 2, 3 and the displacement vector is given by

u = ue1 + ve2 + we3 . (6.23)

The deformation gradient (2.2) in Cartesian coordinates is given by

∂xi
F = Gradx = ei ⊗ Ej i, j = 1, 2, 3. (6.24)
∂Xj
CHAPTER 6. THREE-DIMENSIONAL LARGE DEFORMATIONS 98

alternatively, which is

∂x1 ∂x1 ∂x1

F= e1 ⊗ E1 + e1 ⊗ E2 + e1 ⊗ E3 +
∂X1 ∂X2 ∂X3
∂x2 ∂x2 ∂x2
e2 ⊗ E1 + e2 ⊗ E2 + e2 ⊗ E3 +
∂X1 ∂X2 ∂X3
∂x3 ∂x3 ∂x3
e3 ⊗ E1 + e3 ⊗ E2 + e3 ⊗ E3 . (6.25)
∂X1 ∂X2 ∂X3

6.2.2 Nominal stress and Cauchy stress

The inverse of the deformation gradient F is given by

 

 ∂x2 ∂x3 ∂x2 ∂x3 ∂x1 ∂x3 ∂x1 ∂x3 ∂x1 ∂x2 ∂x1 ∂x2 

 ∂X2 ∂X3 − ∂X3 ∂X2 ∂X3 ∂X2 − ∂X2 ∂X3
 ∂X2 ∂X3 − ∂X3 ∂X2 


 

F−1 = ∂x2 ∂x3 ∂x2 ∂x3
∂X3 ∂X1 − ∂X1 ∂X3
∂x1 ∂x3 ∂x1 ∂x3
∂X1 ∂X3 − ∂X3 ∂X1
∂x1 ∂x2
− ∂x1 ∂x2 (6.26)

 ∂X3 ∂X1 ∂X1 ∂X3 


 


 ∂x2 ∂x3 − ∂x2 ∂x3 ∂x1 ∂x3 − ∂x1 ∂x3 ∂x1 ∂x2 ∂x1 ∂x2 

∂X1 ∂X2 ∂X2 ∂X1 ∂X2 ∂X1 ∂X1 ∂X2 ∂X1 ∂X2 − ∂X2 ∂X1

We now substitute the strain-energy function for neo-Hookean material (6.15) and
(6.26) into the the nominal stress (6.16) for incompressible elastic material to obtain
   

 ∂x ∂x ∂x   ∂x2 ∂x3 ∂x2 ∂x3 ∂x1 ∂x3 ∂x1 ∂x3 ∂x1 ∂x2 ∂x1 ∂x2 
 ∂X1 ∂X1 ∂X1   ∂X2 ∂X3 − ∂X3 ∂X2 ∂X3 ∂X2 − ∂X2 ∂X3 ∂X2 ∂X3 − ∂X3 ∂X2

1 2 3
 
 
 


 
 
 

S = µ ∂X ∂x1 ∂x2 ∂x3 − p ∂x2 ∂x3 ∂x2 ∂x3 ∂x1 ∂x3 ∂x1 ∂x3 ∂x1 ∂x2 ∂x1 ∂x2 .
 2 ∂X2 ∂X2   ∂X3 ∂X1 − ∂X1 ∂X3 ∂X1 ∂X3 − ∂X3 ∂X1 ∂X3 ∂X1 − ∂X1 ∂X3 

 
 
 


   
 ∂x1 ∂x2 ∂x3   
 ∂x2 ∂x3 − ∂x2 ∂x3 ∂x1 ∂x3 − ∂x1 ∂x3 ∂x1 ∂x2 − ∂x1 ∂x2 

∂X3 ∂X3 ∂X3 ∂X1 ∂X2 ∂X2 ∂X1 ∂X2 ∂X1 ∂X1 ∂X2 ∂X1 ∂X2 ∂X2 ∂X1

6.2.3 Governing equations

Using (2.18) we have the equilibrium equations

∂Sij
DivS = Ek · Ei ⊗ ej
∂Xk
∂Sij
= ej = 0 i, j = 1, 2, 3. (6.27)
∂Xi

The final equilibrium equations in component form could be written as

∂S11 ∂S21 ∂S31

+ + = 0, (6.28)
∂X1 ∂X2 ∂X3
∂S12 ∂S22 ∂S32
+ + = 0, (6.29)
∂X1 ∂X2 ∂X3
∂S13 ∂S23 ∂S33
+ + = 0. (6.30)
∂X1 ∂X2 ∂X3

The pressure boundary conditions are given as


 −P F−T N on R = B
ST N = (6.31)

0 on R = A,
CHAPTER 6. THREE-DIMENSIONAL LARGE DEFORMATIONS 99

In component form which are



 −1
S N + S21 N2 = −P (F11 −1
N1 + F21 N2 ),

 11 1
−1 −1 (6.32)
S12 N1 + S22 N2 = −P (F12 N1 + F22 N2 ),



 S N + S N = −P (F −1 N + F −1 N ).
13 1 23 2 13 1 23 2

The end conditions are zero displacements i.e.

u = v = w = 0. (6.33)

And the incompressible condition is

J = det F = 1 (6.34)

and J is

∂x1 ∂x2 ∂x3 ∂x2 ∂x3 ∂x1 ∂x2 ∂x3 ∂x2 ∂x3
J = ( − )+ ( − )
∂X1 ∂X2 ∂X3 ∂X3 ∂X2 ∂X2 ∂X3 ∂X1 ∂X1 ∂X3
∂x1 ∂x2 ∂x3 ∂x2 ∂x3
+ ( − ). (6.35)
∂X3 ∂X1 ∂X2 ∂X2 ∂X1

6.2.4 Unit outward normal to the outer surface

Let P be the arbitrary point on the outer surface, O be the original point located on the
center of the bottom of the cylinder. We choose the outer radius of the cylinder B = 1.
−−→
We use R denoting the position vector OP which is given by

R = B cos ΘE1 + B sin ΘE2 + ZE3 (6.36)

The unit outward normal to the outer surface is given by

RΘ × RZ
N=
k RΘ × RZ k
= cos ΘE1 + sin ΘE2
= X1 E1 + X2 E2 . (6.37)

6.3 Linear case

The linear case is presented for comparison with the above nonlinear one.
CHAPTER 6. THREE-DIMENSIONAL LARGE DEFORMATIONS 100

6.3.1 Deformation gradient and Cauchy stress

The deformation gradient in matrix form,

 

 ∂u1 ∂u1 ∂u1 


 ∂x1 + 1 ∂x2 ∂x3 


 

F= ∂u2 ∂u2 ∂u2 . (6.38)
+1

 ∂x1 ∂x2 
∂x3


 


 ∂u3 
∂x1
∂u3
∂x2
∂u3
∂x3 +1

We have the Cauchy stress tensor for isotropic incompressible material in the form

σ = −pI + µ[gradu + (gradu)T ], (6.39)

where

∂ui
gradu = ei ⊗ ej i, j = 1, 2, 3. (6.40)
∂xj

The final component form of Cauchy stress is given by

 

 ∂u p ∂u ∂u ∂u ∂u 
 2 ∂x1 − µ ∂x2 + ∂x1 ∂x3 + ∂x1 
1 1 2 1 3
 


 

∂u ∂u ∂u p ∂u
σ = µ ∂x1 + ∂x2 2 ∂x2 − µ ∂x2 + ∂x3 ∂u (6.41)

 2 1 2 3 2 


 


 ∂u1 + ∂u3 ∂u2 + ∂u3 2 ∂u3 − p  
∂x3 ∂x1 ∂x3 ∂x2 ∂x3 µ

6.3.2 Governing equations

Using the equilibrium equation (2.19) we obtain

∂σ11 ∂σ21 ∂σ31

+ + = 0, (6.42)
∂x1 ∂x2 ∂x3
∂σ12 ∂σ22 ∂σ32
+ + = 0, (6.43)
∂x1 ∂x2 ∂x3
∂σ13 ∂σ23 ∂σ33
+ + = 0. (6.44)
∂x1 ∂x2 ∂x3

And the incompressible condition is

∂ui ∂u1 ∂u2 ∂u3

divu = = + + = 0, (6.45)
∂xi ∂x1 ∂x2 ∂x3

where u is the displacement vector.

The pressure boundary conditions are

 −P n on R = B
σn = (6.46)

0 on R = A,
CHAPTER 6. THREE-DIMENSIONAL LARGE DEFORMATIONS 101

where n the unit outer normal to the external lateral surface of the tube. The component
form is given by


 σ n + σ12 n2 = −P n1 ,

 11 1
σ21 n1 + σ22 n2 = −P n2 , (6.47)



 σ n + σ n = 0.
31 1 32 2

The end conditions are zero displacement i.e.

u = v = w = 0. (6.48)
Chapter 7

Conclusions

The main contributions of this research are the investigation of both the axisymmetric
and asymmetric bifurcations of thick-walled circular cylindrical elastic tubes under axial
loading and external pressure based on the work of Haughton and Ogden [35] and the
nonlinear analysis of the finite axisymmetric deformations of an elastic tube under external
pressure based on the general fully nonlinear governing equations in Lagrangian form, and
finally, the development of the partial differential equations for the three dimensional large
deformations of an elastic tube under external pressure.
As stated in the discussion and conclusion in Chapters 4 and 5, we summarize our
main results here.
In the work [78] (in Chapter 4), we found that the mode-2 bifurcation pressure is
insensitive to the variation of the wall thickness and deviate from the thin shell prediction
in both the very thin and thick-walled regimes. The mode transition from lower to higher
ones occurs in the range of axial compression, for very short and sufficiently thick tubes.
We have also shown that, contrary to thin-shell theory, for sufficiently thick tubes, mode
2 bifurcation becomes the only dominant mode, without transition from lower to higher
modes for sufficiently short tubes.
In the work [79] (in Chapter 5), as expected, the linear and nonlinear models represent
nearly the same results for small deformation. However, large differences have been found
in the predictions of the linear and nonlinear models under large external pressure, and
the dominant nonlinear features are the corner bulging, for relatively short tubes and,
the occurrence of higher mode deformations for longer tubes. It has been observed that
material nonlinearity dominates the deformations in the shorter and thicker tubes, for
which the strains computed are larger, while geometrical nonlinearity seems to provide

102
CHAPTER 7. CONCLUSIONS 103

more influence on the deformations in the longer and thinner tubes, for which the strains
are much smaller. The shear components of Cauchy stress exhibits the greatest differences
between the results of the linear and nonlinear theories. Shear splitting, with variation
in signs of the shear stress in neighbouring regions is a unique nonlinear feature. The
boundary layer formation at the ends of the long, thin tubes, which is a hallmark of classical
shell theory, is also represented in the nonlinear analysis. This is the first systematic
nonlinear analysis of externally pressurized thick-walled elastic tubes, although the simple
neo-Hookean material has been adopted, and the results may have significant implications
for certain physiological problems involving soft vessels undergoing large deformation.
Regarding future projects, there are several possibilities as follows.
In the next phase of the work in Chapter 4 we shall investigate the post-buckling
behaviour of elastic tubes under external pressure and axial loading. In particular, the
effect of wall-thickness on compliance of the tubes between buckling and self contact will
be studied in order to interpret the puzzling phenomenon that for tubes subjected to
external pressure, after a certain degree of collapse, thick tubes may be more compliant
than thinner ones [11, 52].
Since we have dealt with axisymmetric problems in Chapter 5, we can only simu-
late the necked or barrelled states of a cylinder. In addition, we have not carried out
any bifurcation analysis and it remains to ascertain the stability status of the solutions
obtained, although previous studies on similar nonlinear problems, albeit with different
boundary conditions [55], have shown that there exist nontrivial axisymmetric stable (half
neck or multiple-neck) solutions. However, in many physiological situations, nonlinear
deformations that break this symmetry (both for the original deformation or the bifur-
cation analysis) could be more significant. We shall continue to pursue this problem in
subsequent work based on the results of Chapter 5. Considering the fully nonlinear struc-
tures of the system of the equations governing the three dimensional deformations, in the
numerical solution, resulting in a dense and nonsymmetric coefficient matrix it is still not
sure that the numerical simulation could be carried out successfully using the present finite
element method and the algorithms for solving nonlinear sets of equations. (see details
in section 9 in the book [62]). In order to avoid the difficulty of solving the nonlinear
systems directly, we note that a more natural and effective analysis approach has been
developed by Bathe [8] by referring all variables to a previously known calculated equi-
librium configuration and linearizing the resulting equations to obtain an approximate
CHAPTER 7. CONCLUSIONS 104

solution. Once the above questions are answered, a natural next step would be to replace
the zero-displacement end boundary condition by axial loading including both compres-
sion and extension to provide a better approximation of the real physiological problem as
mentioned in the Introduction.
Another very interesting but difficult future line of research would be to extend the
formulations of all the above problems in Chapter 4-6 to the corresponding dynamics cases.
Appendix A

Appendix A

In this appendix we represent the derivation of discretizing integrations for the nonlinear
axisymmetric case in Chapter 5. Note that the subscripts of components of nominal stress
tensor S are changed. The connection between the alternative notations and the one used
in Chapter 5 are S11 = SRr , S31 = SZr , S13 = SRz , S33 = SZz , S22 = SΘθ .

A.1 Integration of the equilibrium equations by parts

Applying Galerkin’s method to the differential equilibrium equations,

Z
1
Ni [S11,1 + S31,3 + (S11 − S22 )]dΩ = 0 (A.1)
Ω Z R
1
Ni (S13,1 + S33,3 + S13 )dΩ = 0 (A.2)
Ω R
where Ni is the i th shape function and Ω is the volume of the cylinder in reference
configuration.
For axisymmetric deformations, the nominal stress components and shape functions
are independent of the coordinate Θ, we can decompose the integral above as follows,
Z
1
RNi [S11,1 + S31,3 + (S11 − S22 )]dRdZ = 0 (A.3)
R,Z R
Z
1
RNi (S13,1 + S33,3 + S13 )dRdZ = 0 (A.4)
R,Z R
Integration by parts, the first term of (A.3) can be written as
Z Z Z
RNi S11,1 dRdZ = − (RNi ),1 S11 dRdZ + (RNi S11 )|R
R1 dZ
2

R,Z R,Z Z

The second term of (A.3) can be written as

Z Z Z
RNi S31,3 dRdZ = − RNi,3 S31 dRdZ + R(Ni S31 )|Z
Z1 dR
2

R,Z R,Z R

105
APPENDIX A. APPENDIX A 106

and hence we obtain from (A.3),

Z Z Z
R2
[R(Ni,1 S11 + Ni,3 S31 ) + Ni S22 ]dRdZ = (RNi S11 )|R1 dZ + R(Ni S31 )|Z
Z1 dR
2
(A.5)
R,Z Z R

The similar techniques operate on (A.4), we have

Z Z Z
R(Ni,1 S13 + Ni,3 S33 )dRdZ = (RNi S13 )|R
R1
2
dZ + R(Ni S33 )|Z
Z1 dR (A.6)
2

R,Z Z R

We use notations int1 and int2 to denote the two integrands in (A.5) and (A.6),
respectively,

int1 = R(Ni,1 S11 + Ni,3 S31 ) + Ni S22

int2 = R(Ni,1 S13 + Ni,3 S33 )

By the help of mathematica, we obtain

int1 = t2 Ni
½· ¸ ¾
(1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
+ RNi,1 + (1 + uR ) + sin 2ψ( − )uZ
½ 2λ 1 2λ 3 2 λ1 λ3
1 t1 t3
+ RNi,3 sin 2ψ( − )(1 + uR )
· 2 λ 1 λ3 ¸ ¾
(1 − cos 2ψ)t1 (1 + cos 2ψ)t3
+ + uZ (A.7)
2λ1 2λ3

First term of (A.7),

Z Z
t2 Ni dRdZ = Ni (µλ2 − pλ−1
2 )dRdZ
R,Z ZR,Z
u R
= Ni (µ(1 + ) − p)dRdZ
R R+u
ZR,Z Z Z
µ R 0
= µNi dRdZ + Ni Nj dRdZuj − Ni Nj dRdZpj
R,Z R,Z R R,Z R + u

Using = µ − pλ−2
t1
λ1 1 , second term of (A.7),
Z ½· ¸ ¾
(1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
RNi,1 + (1 + uR ) + sin 2ψ( − )uZ dRdZ
R,Z 2λ1 2λ3 2 λ1 λ3
Z · ¸
(1 + cos 2ψ) (1 − cos 2ψ)t 3
= RNi,1 (µ − λ−21 p) + dRdZ
R,Z 2 2λ 3
Z · ¸ Z
(1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
+ RNi,1 + uR dRdZ + RNi,1 sin 2ψ( − )uZ dRdZ
2λ1 2λ 3 2 λ1 λ3
ZR,Z · ¸ R,Z
1 + cos 2ψ (1 − cos 2ψ)t3
= RNi,1 µ+ dRdZ
2 2λ3
Z R,Z ½· ¸ ¾
(1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
+ RNi,1 + Nj,1 + sin 2ψ( − )Nj,3 dRdZuj
R,Z 2λ1 2λ3 2 λ1 λ3
Z
1 + cos 2ψ −2 0
− R λ1 Ni,1 Nj dRdZpj
R,Z 2
APPENDIX A. APPENDIX A 107

Third term of (A.7),

Z ½ · ¸ ¾
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
RNi,3 sin 2ψ( − )(1 + uR ) + + uZ dRdZ
R,Z 2 λ1 λ3 2λ1 2λ3
Z
1 t1 t3
= RNi,3 sin 2ψ( − )dRdZ
2 λ1 λ3
Z R,Z ½ · ¸ ¾
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
+ RNi,3 sin 2ψ( − )Nj,1 + + Nj,3 dRdZuj
R,Z 2 λ1 λ3 2λ1 2λ3
Hence from equation (A.5), we obtain
Z ½· ¸ ¾
µ (1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
Ni Nj + RNi,1 + Nj,1 + sin 2ψ( − )Nj,3
R,Z R ½ 2λ1
·
2λ3 2
¸
λ
¾1
λ3
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
+ RNi,3 sin 2ψ( − )Nj,1 + + Nj,3 dRdZuj
Z ·2 λ1 λ3 ¸2λ1 2λ
Z 3
1 1 + cos 2ψ −2 0
− R Ni + λ1 Ni,1 Nj dRdZpj = − µNi dRdZ (A.8)
R+u 2
ZR,Z · ¸ Z R,Z
1 + cos 2ψ (1 − cos 2ψ)t3 1 t1 t3
− RNi,1 µ+ dRdZ − RNi,3 sin 2ψ( − )dRdZ
2 2λ3 R,Z 2 λ1 λ3
ZR,Z Z
+ (RNi S11 )|R R1 dZ +
2
R(Ni S31 )|Z
Z1 dR
2

Z R

And the integrand of (A.6),

½· ¸ ¾
(1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
int2 = RNi,1 + wR + sin 2ψ( − )(1 + wZ )
½ 2λ1 2λ3· 2 λ1 λ3 ¸ ¾
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
+ RNi,3 sin 2ψ( − )wR + + (1 + wZ(A.9)
)
2 λ1 λ3 2λ1 2λ3
First term of (A.9),
Z ½· ¸ ¾
(1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
RNi,1 + wR + sin 2ψ( − )(1 + wZ ) dRdZ
R,Z 2λ1 2λ3 2 λ1 λ3
Z
1 t1 t3
= RNi,1 sin 2ψ( − )dRdZ
2 λ1 λ3
Z R,Z ½· ¸ ¾
(1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
+ RNi,1 + Nj,1 + sin 2ψ( − )Nj,3 dRdZwj
R,Z 2λ1 2λ3 2 λ1 λ3

Using = µ − pλ−2
t3
λ3 3 , second term of (A.9),
Z ½ · ¸ ¾
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
RNi,3 sin 2ψ( − )wR + + (1 + wZ ) dRdZ
R,Z 2 λ1 λ3 2λ1 2λ3
Z · ¸
(1 − cos 2ψ)t1 1 + cos 2ψ
= RNi,3 + (µ − λ−2
3 p) dRdZ
R,Z 2λ 1 2
Z ½ · ¸ ¾
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
+ RNi,3 sin 2ψ( − )Nj,1 + + Nj,3 dRdZwj
2 λ1 λ3 2λ1 2λ3
ZR,Z · ¸
(1 − cos 2ψ)t1 1 + cos 2ψ
= RNi,3 + µ dRdZ
2λ1 2
Z R,Z ½ · ¸ ¾
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
+ RNi,3 sin 2ψ( − )Nj,1 + + Nj,3 dRdZwj
R,Z 2 λ1 λ3 2λ1 2λ3
Z
(1 + cos 2ψ) −2 0
− R λ3 Nj,3 Nj dRdZpj
R,Z 2
APPENDIX A. APPENDIX A 108

Hence from (A.6), we obtain

Z ½· ¸ ¾
(1 + cos 2ψ)t1 (1 − cos 2ψ)t3 1 t1 t3
RNi,1 + Nj,1 + sin 2ψ( − )Nj,3
2λ1 2λ 2 λ1 λ3
R,Z ½ · 3 ¸ ¾
1 t1 t3 (1 − cos 2ψ)t1 (1 + cos 2ψ)t3
+ RNi,3 sin 2ψ( − )Nj,1 + + Nj,3 dRdZwj
Z 2 λ1 λ3 2λZ1 2λ3
(A.10)
1 + cos 2ψ −2 0 1 t1 t3
− R λ3 Ni,3 Nj dRdZpj = − RNi,1 sin 2ψ( − )dRdZ
2 2 λ1 λ3
ZR,Z · ¸R,Z Z Z
(1 − cos 2ψ)t1 1 + cos 2ψ
− RNi,3 + µ dRdZ + (RNi S13 )|R 2
R1 dZ + R(Ni S33 )|Z
Z1 dR
2

R,Z 2λ1 2 Z R

Incompressible condition

λ1 λ3 = λ−1
2

expanding incompressible condition above, we have

1
u + uR + (1 + uR )wZ − uZ wR = 0 (A.11)
R+u

Integration (A.11), we obtain

Z · ¸
0 1
RNi u + uR + (1 + uR )wZ − uZ wR dRdZ
R,Z R+u
Z µ ¶ Z
0 1 0
= RNi Nj + Nj,1 dRdZuj + RNi [(1 + uR )Nj,3 − uZ Nj,1 ] dRdZwj
R,Z R+u R,Z

Finally, we have
Z µ ¶
0 1
RNi Nj + Nj,1 dRdZuj (A.12)
R,Z R+u
Z
0
+ RNi [(1 + uR )Nj,3 − uZ Nj,1 ] dRdZwj = 0
R,Z
Appendix B

Appendix B

In this appendix we present the mesh files used in Chapter 5. For details of the format of
these files, we refer to the web page of libmesh: http://libmesh.sourceforge.net/publications.php.

B.1 Mesh files

The file mesh.xda for a tube with A/B = 0.5, L/B = 1 is give by
LIBM 0
4 # Num. of elements
6 # Num. nodes
20 # length of connectivity vector
6 # Num. boundary conds
65536 # string size
1 # Num. elements blocks
3 # Element types in each block
4 # Num. of elements in each block at each level

Id String
Title String
0 1 2 0 -1
0 2 3 1 -1
3 2 5 2 -1
2 4 5 3 -1

109
APPENDIX B. APPENDIX B 110

0.5 0 0
100
1 0.5 0
0.5 0.5 0
110
0.5 1 0

000
011
122
222
301
313

The file mesh.xda for a tube with A/B = 0.5, L/B = 5 is give by
LIBM 0
20 # Num. of elements
22 # Num. nodes
100 # length of connectivity vector
22 # Num. boundary conds
65536 # string size
1 # Num. elements blocks
3 # Element types in each block
20 # Num. of elements in each block at each level

Id String
Title String
0 1 2 0 -1
0 2 3 1 -1
3 2 4 2 -1
3 4 5 3 -1
5 4 6 4 -1
5 6 7 5 -1
APPENDIX B. APPENDIX B 111

7 6 8 6 -1
7 8 9 7 -1
9 8 10 8 -1
9 10 11 9 -1
11 10 13 10 -1
10 12 13 11 -1
13 12 15 12 -1
12 14 15 13 -1
15 14 17 14 -1
14 16 17 15 -1
17 16 19 16 -1
16 18 19 17 -1
19 18 21 18 -1
18 20 21 19 -1

0.5 0 0
100
1 0.5 0
0.5 0.5 0
110
0.5 1 0
1 1.5 0
0.5 1.5 0
120
0.5 2 0
1 2.5 0
0.5 2.5 0
130
0.5 3 0
1 3.5 0
0.5 3.5 0
140
0.5 4 0
APPENDIX B. APPENDIX B 112

1 4.5 0
0.5 4.5 0
150
0.5 5 0

000
011
122
211
322
411
522
611
722
811
922
10 2 2
11 0 1
12 2 2
13 0 1
14 2 2
15 0 1
16 2 2
17 0 1
18 2 2
19 0 1
19 1 3
References

[1] M. Amabili and M.P. Paı̈doussis. Review of studies on geometrically nonlinear vi-
brations and dynamics of circular cylindrical shells and panels, with and without
fluid-structure intersection. Applied Mechanics Reviews, 56:349–381, 2003.

[2] M. Anliker, W.E. Moritz, and E. Ogden. Transmission characteristics of axial waves
in blood vessels. Journal of Biomechanics, 1:235–246, 1968.

[3] B. Bader, R. Hooper, T. Kolda, R. Pawlowski, E. Phipps, and A. Salinger. Nox

and loca: Object-oriented nonlinear solver and continuation packages. The Trilinos
Project website http://trilinos.sandia.gov/.

[4] S. Balay, K. Buschelman, V. Eijkhout, W.D. Gropp, D. Kaushik, M.G. Knepley,

L.C. McInnes, B.F. Smith, and H. Zhang. Petsc users manual. Technical Report
ANL-95/11 – Revision 3.0.0, Argonne National Laboratory, 2008.

[5] F.C. Bardi and S. Kyriakides. Plastic buckling of circular tubes under axial
compression–part i: Experiments. International Journal of Mechanical Sciences,
48:830–841, 2006.

[6] F.C. Bardi, S. Kyriakides, and H.D. Yun. Plastic buckling of circular tubes under
axial compression–part ii: Analysis. International Journal of Mechanical Sciences,
48:842–854, 2006.

[7] S.B. Batdorf. A simplified method of elastic stability analysis for thin cylindrical
shells. NACA Report, page 874, 1947.

[8] K.J. Bathe. Finite ELement Procedures. Prentice Hall, New Jersey, 1996.

[9] M.F. Beatty and P. Dadras. Some experiments on the elastic stability of some highly
elastic bodies. International Journal of Engineering Science, 14:233–238, 1968.

113
REFERENCES 114

[10] C.D. Bertram. Two modes of instability in a thick-walled collapsible tube conveying
a flow. Journal of Biomechanics, 15:223–224, 1982.

[11] C.D. Bertram. The effects of wall thickness, axial strain and end proximity on the
pressure-area relation of collapsible tubes. Journal of Biomechanics, 20:863–876, 1987.

[12] C.D. Bertram and N.S.J. Elliot. Flow-rate limitation in a uniform thin-walled col-
lapsible tube, with comparison to a uniform thick-walled tube and a tube of tapering
thickness. Journal of Fluids and Structures, 17:541–559, 2003.

[13] C.D. Bertram, C.J. Raymond, and T.J. Pedley. Mapping of instabilities for flow
through collapsed tubes of differing length. Journal of Fluids and Structures, 4:125–
153, 1990.

[14] C.D. Bertram, C.J. Raymond, and T.J. Pedley. Application of dynamical system
concepts to the analysis of self-excited oscillations of a collapsible tube conveying a
flow. Journal of Fluids and Structures, 5:391–426, 1991.

[15] G. Buzzi-Ferraris. Scientific C++: Building Numerical Libraries the Object-Oriented

Way. Addison-Wesley, Reading, MA, 1993.

[16] G. Buzzi-Ferraris. Object-Oriented Analysis and Design with Applications. Benjamin-

Cummings, Redwood City, CA, 1994.

[17] Z.X. Cai and X.Y. Luo. A fluidbeam model for flow in collapsible channel. Journal
of Fluids and Structures, 17 (1):123–144, 2003.

[18] D.M. Cannell. George Green mathematician and physicist. The Athlone Press, Lon-
don, 1993.

[19] G.F. Carey, M. Anderson, B. Carnes, and B. Kirk. Some aspects of adaptive grid tech-
nology related to boundary and interior layers. J. Comput. Appl. Math., 166(1):5586,
2004.

[20] H. Demiray. Nonlinear waves in a prestressed thick elastic tube filled with an invisid
fluid. International Journal of Engineering Science, 34:1519–1530, 1996.

[21] H. Demiray. Nonlinear wave modulation in a fluid filled thick elastic tube. Interna-
tional Journal of Engineering Science, 36:1061–1082, 1998.
REFERENCES 115

[22] Hinton E. and Campbell J. Local and global smoothing of discontinuous finite element
function using a least squares method. Int. J. Numer. Meth. Eng., 8:461–480, 1974.

[23] H.A. Erbay and H. Demiray. Finite axisymmetric deformations of elastic tubes: An
approximate method. Journal of Engineering Mathematics, 29:451–472, 1995.

[24] J.E. Flaherty, J.B. Keller, and S.I. Rubinow. Post buckling behavior of elastic tubes
and rings with opposite sides in contact. SIAM Journal of Applied Mathematics,
23:446–455, 1972.

[25] W. Flügge. Stresses in Shells. Springer, Berlin, second edition, 1973.

[26] Y.B. Fu and R.W. Ogden. Nonlinear stability analysis of pre-stressed elastic bodies.
Continuum Mech. Thermodyn., 11:141172, 1999.

[27] Y.B. Fu and R.W. Ogden. Nonlinear Elasticity: Theroy and Application, volume 283
of London Mathematical Society Lecture Notes Series. Cambridge University Press,
2001.

[28] A.E. Green, R.S. Rivilin, and R.T. Shield. General theory of small elastic deformations
superposed on finite elastic deformations. Proceedings of the Royal Society of London
A, 211:128–154, 1952.

[29] D.E. Gregg and L.C. Fisher. Blood supply to the heart. In Handbook of Physiology
(ed. Hamilton W.F. and Dow P.). American Physiological Society, Washington D.C.,
1963.

[30] A.C. Guyton and L.H. Adkins. Quantitative aspects of the collapse factor in relation
to venous return. American Journal of Physiology, 177:523–527, 1954.

[31] Si H. A quality tetrahedral mesh generator and a 3d delaunay triangulator.

http://tetgen.berlios.de/.

[32] D.M. Haughton and R.W. Ogden. On the incremental equations in non-linear
elasticity–i. membrane theory. Journal of the Mechanics and Physics of Solids, 26:93–
110, 1978a.

[33] D.M. Haughton and R.W. Ogden. On the incremental equations in non-linear
elasticity–ii. bifurcation of pressured spherical shells. Journal of the Mechanics and
Physics of Solids, 26:111–138, 1978b.
REFERENCES 116

[34] D.M. Haughton and R.W. Ogden. Bifurcation of inflated circular cylinders of elastic
material under axial loading–i. membrane theory for thin-walled tubes. Journal of
the Mechanics and Physics of Solids, 27:179–212, 1979a.

[35] D.M. Haughton and R.W. Ogden. Bifurcation of inflated circular cylinders of elastic
material under axial loading–ii. exact theory for thick-walled tubes. Journal of the
Mechanics and Physics of Solids, 27:489–512, 1979b.

[36] M. Heil. The stability of cylindrical shells conveying viscous flow. Journal of Fluids
and Structures, 10:173–196, 1996.

[37] M. Heil. Stokes flow in collapsible tubes: computation and experiment. Journal of
Fluid Mechanics, 353:258–312, 1997.

[38] M. Heil and T.J. Pedley. Large axisymmetric deformation of a cylindrical shell con-
veying a viscous flow. Journal of Fluids and Structures, 9:237–256, 1995.

[39] M. Heil and T.J. Pedley. Large post-buckling deformations of cylindrical shells con-
veying viscous flow. Journal of Fluids and Structures, 10:565–599, 1996.

[40] B.P.C. Ho and S. Cheng. Some problems in stability of heterogeneous aeolotropic

cylindrical shells under combined loading. AIAA Journal, 1:1603–1607, 1963.

[41] O.E. Jensen and M. Heil. High-frequency self-excited oscillations in a collapsible-

channel flow. Journal of Fluid Mechanics, 481:235–268, 2003.

[42] R.D. Kamm and T.J. Pedley. Flow in collapsible tubes: a brief review. ASME Journal
of Biomechanical Engineering, 111:177–179, 1989.

[43] D.W. Kelly, J. P. Gago, O.C. Zienkiewicz, and I. Babuska. A posteriori error analysis
and adaptive processes in the finite element method: Part i error analysis. Int. J.
Num. Meth. Engng., 19:15931619, 1983.

[44] B. Kirk, J.W. Peterson, R.H. Stogner, and G.F. Carey. Libmesh: A c++ library for
parallel adaptive mesh refinement/coarsening simulations. Engineering with Comput-
ers, 22:237–254, 2006.

[45] B.S. Kirk. Adaptive Finite Element Simulation of Flow and Transport Applications
on Parallel Computers. The University of Texas at Austin, 2007.
REFERENCES 117

[46] A.N. Kounadis. Recent advances on postbuckling analyses of thin-walled struc-

tures: Beams, frames and cylindrical shells. Journal of Constructional Steel Research,
62:1101–1115, 2006.

[47] A. Libai and J.G. Simmonds. The Nonlinear Theory of Elastic Shells. Cambridge
University Press, Cambridge, 1998.

[48] X.Y. Luo and T.J. Pedley. A numerical simulation of steady flow in a 2d collapsible
channel. Journal of Fluids and Structures, 9:149–174, 1995.

[49] X.Y. Luo and T.J. Pedley. A numerical simulation of unsteady flow in a 2d collapsible
channel. Journal of Fluid Mechanics, 314:191–225, 1996.

[50] X.Y. Luo and T.J. Pedley. The effects of the wall inertia on the 2d collapsible channel
flow. Journal of Fluid Mechanics, 363:253–280, 1998.

[51] X.Y. Luo and T.J. Pedley. Flow limitation and multiple solutions in 2d collapsible
channel flow. Journal of Fluid Mechanics, 420:301–324, 2000.

[52] A. Marzo, X.Y. Luo, and C.D. Bertram. Three-dimensional collapse and steady flow
in thick-walled flexible tubes. Journal of Fluids and Structures, 20:817–835, 2005.

[53] T.B. Moodie and G.E. Swaters. Nolinear waves and shock calculations for a hypere-
lastic fluid-filled tube. Quarterly Applied Mathematics, 47:705–732, 1989.

[54] W.A. Nash. Buckling of thin cylindrical shells subject to hydrostatic pressure. Journal
of the Aeronautical Sciences, 21:354–355, 1954.

[55] P.V. Negrón-Marrero. An analysis of the linearized equations for axisymmetric defor-
mations of hyperelastic cylinders. Mathematics and Mechanics of Solids, 4:109–133,
1999.

[56] J.L. Nowinski and M. Shahinpoor. Stability of an elastic circular tube of arbitrary
wall thickness subjected to an external pressure. International Journal of Non-Linear
Mechanics, 4:143–158, 1969.

[57] R.W. Ogden. Large deformation isotropic elasticity i: On the correlation of theory
and experiment for incompressible rubberlike solids. Proceedings of the Royal Society
of London A, 326:565–584, 1972.
REFERENCES 118

[58] R.W. Ogden. On stress rates in solid mechanics with application to elasticity theory.
Proceedings of the Cambridge Philosophical Society, 75:303–319, 1974.

[59] R.W. Ogden. Elastic deformation of rubberlike solids. Mechanics of Solids, The
Rodney Hill 6oth Anniversary Volume:499–537, 1982.

[60] R.W. Ogden. Non-linear Elastic Deformations. Dover Publications, New York, 1997.

[61] T.J. Pedley. Longitudinal tension variation in collapsible channels: a new mechanism
for the breakdown of steady flow. ASME Journal of Biomechanical Engineering,
114:60–67, 1992.

[62] W.H. Press, S.A. Teukolsky, W.T. Vetterling, and B.P. Flannery. Numerical Recipes
in Fortran The Art of Scientific Computing Second Edition. Cambridge University
Press, Cambridge, 1992.

[63] S. Rodbard and L.H. Adkins. A hydrodynamics mechanism for autoregulation of flow.
Cardiologia, 48:532–535, 1966.

[64] G. Rudinger. Shock waves in a mathematical model of the aorta. Journal of Applied
Mechanics, 37:34–37, 1970.

[65] R.A. Scroggs, E.A. Patterson, and S.B.M. Beck. Numerical fluidstructure interaction
model of a collapsible tube using ls-dyna. IUTAM Symposium on Flowin Collapsible
Tubes and Past other Highly Compliant Boundaries, Warwick, UK.

[66] L.H. Sobel. Effects of boundary conditions on the stability of cylinders subject to
lateral and axial pressure. AIAA Journal, 2:1437–1440, 1964.

[67] R.H. Stogner. Parallel Adaptive C 1 Macro-Elements for Nonlinear Thin Film and
Non-Newtonian Flow Problems. The University of Texas at Austin, 2008.

[68] B. Stroustrup. The C++ Progamming Language. Addison-Wesley, 2007.

[69] R. J. Tait, D. J. Steigmann, and J. L. Zhong. Finite twist and extension of a cylindrical
elastic membrane. Acta Mechanica, 117(1):129–143, 1996.

[70] D. Tang and C. Yang. 3-d steady and unsteady blood flowin stenotic collapsible
arteries with symmetric and asymmetric stenoses. Recent Advances in Biomechan-
ics,Proceedings of First International Young Investigators Workshop on Biomechan-
ics, Beijing, 121:171191, 2001.
REFERENCES 119

[71] D. Tang, J. Yang, and D.N. Ku. A nonlinear axisymmetric model with fluidwall in-
teractions for viscous flows in stenotic elastic tubes. ASME Journal of Biomechanical
Engineering, 121:494501, 1999.

[72] R. von Mises. Der kritische außendruck zylindrischer rohre. VDI Zeitschrift, 58:750–
755, 1914.

[73] A.S.D. Wang and A. Ertepinar. Stability and vibrations of elastic thick-walled cylin-
drical and spherical shells subjected to pressure. International Journal of Non-Linear
Mechanics, 7:539–555, 1972.

[74] M. Weissman and L. Mockros. The mechanics of a collapsing tube heart pump.
International Journal of Mechanical Sciences, 9:113–121, 1967.

[75] N. Yamaki. Buckling of circular cylindrical shells under external pressure. Reports of
the Institute of High Speed Mechanics, 20:35–55, 1969.

[76] N. Yamaki. Influence of prebuckling deformation on the buckling of circular cylindrical

shells under external pressure. Reports of the Institute of High Speed Mechanics,
21:81–104, 1970.

[77] N. Yamaki. Elastic Stability of Circular Cylindrical Shells. North-Holland, Amster-

dam, 1984.

[78] Y. Zhu, X.Y. Luo, and R.W. Ogden. Asymmetric bifurcations of thick-walled circular
cylindrical elastic tubes under axial loading and external pressure. International
Journal of Solids and Structures, 45:3410–3429, 2008.

[79] Y. Zhu, X.Y. Luo, and R.W. Ogden. Nonlinear axisymmetric deformations of an
elastic tube under external pressure. European Journal of Mechanics - A/Solids, In
Press, Accepted Manuscript:doi:10.1016/j.euromechsol.2009.10.004, 2010.

[80] O.C. Zienkiewicz, R.L. Taylor, and J.Z Zhu. The finite element method its basis and
fundamentals. Elsevier Butterworth-Heinemann, Oxford, 2005.