CONTINUUM-BASED FINITE ELEMENTS FOR SHAKEDOWN …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/... · the computation of collapse loads and shakedown limits of pipelines and pressure

CONTINUUM-BASED FINITE ELEMENTS FOR SHAKEDOWN ANALYSIS OF

PIPELINES AND PRESSURE VESSELS

Ricardo Rodrigues Martins

Tese de Doutorado apresentada ao Programa de

Pós-graduação em Engenharia Mecânica,

COPPE, da Universidade Federal do Rio de

Janeiro, como parte dos requisitos necessários à

obtenção do título de Doutor em Engenharia

Mecânica.

Orientadores: Nestor Alberto Zouain Pereira

Lavinia Maria Sanábio Alves Borges

Rio de Janeiro

Setembro de 2013

iii

Martins, Ricardo Rodrigues

Continuum-based finite elements for shakedown

analysis of pipelines and pressure vessels/ Ricardo

Rodrigues Martins. – Rio de Janeiro: UFRJ/COPPE,

2013.

XV, 138 p.: il.; 29,7 cm.


Lavinia Maria Sanábio Alves Borges

Tese (doutorado) – UFRJ/ COPPE/ Programa de

Engenharia Mecânica, 2013.

Referências Bibliográficas: p. 127-138.

1. Direct methods. 2. Shakedown. 3. Finite Element

Method. I. Pereira, Nestor Alberto Zouain, et al. II.

Universidade Federal do Rio de Janeiro, COPPE,

Programa de Engenharia Mecânica. III. Título.

To God and my family

iv

Acknowledgments

I would like to express my sincere gratitude to all who contributed to this work. First

of all, I fully acknowledge the fundamental support given by Petroleo Brasileiro

S.A. (PETROBRAS) which believes in the potential of its collaborators, giving

us excellent opportunities for continuous learning. In particular, I would like to

thank Petrobras Research and Development Center (CENPES ), in the person of

the managers, Arthur Curty Saad and Luiz Augusto Petrus Levy, and all those

from PETROBRAS who trusted in this project from the very beginning and made

it possible the conduction of this research.

I am fully thankful to my advisor, Prof. Nestor Zouain, for his availability,

guidance, patience and friendship. This work can be regarded in great part as a

result of his valuable ideas. All the steps of this thesis were carried out under

his careful supervision, and I am quite sure that this work would not be the same

without his support.

In a similar way, I would like to express my gratitude to my co-advisor, Prof.

Lavinia Borges. She provided excellent contributions in different phases of this work.

Her enthusiasm helped me to be always motivated, especially during the hard times.

I am totally grateful to Prof. Eduardo de Souza Neto, who accepted to be my

supervisor during my stay as visiting researcher at Swansea University. His contri-

butions and insights were very important and improved the final results of the thesis.

Moreover, his friendship created an exceptional environment to the conduction of

part of this work in Wales.

I would like to express my gratitude to my colleagues from Pipelines and Risers

Department (CENPES/PDEP/TDUT), especially those who I have been working

directly: Adilson Benjamin, Ana Lucia Torres, Dary Kayser Jr., Dauro Noronha

Jr., Divino Cunha, Edmundo Queiroz, Erico Santos, Lea Troina, Ludimar Aguiar ,

Marcio Mourelle, Marcos Martins, Stael Senra. I am fully grateful for their support

and encouragement. I am also thankful for the support from the secretaries Arleide

Araujo, Valeria Santos and Andreia Albuquerque.

I would like to show my gratitude to the faculty members and staff from Graduate

Program in Mechanical Engineering, COPPE/PEM, of Federal University of Rio de

Janeiro for sharing their knowledge and time.

v

Finally, but not less importantly, I am very grateful to my wife Daiane, my

parents Wilson and Lisete, and my sisters Priscila and Pamela, for their true love

and unconditional support. They give me reason for trying to be a bit better every

day.

vi

Resumo da Tese apresentada a COPPE/UFRJ como parte dos requisitos necessarios

para a obtencao do grau de Doutor em Ciencias (D.Sc.)

ELEMENTOS FINITOS BASEADOS NO CONTINUO PARA ANALISE DE

ACOMODACAO ELASTICA DE DUTOS E VASOS DE PRESSAO


Setembro/2013


Lavinia Maria Sanabio Alves Borges

Programa: Engenharia Mecanica

Este trabalho trata do desenvolvimento de elementos finitos estruturais baseados

no contınuo para o calculo por metodos diretos do colapso e da acomodacao elastica

em dutos e vasos de pressao. Primeiramente e realizada uma revisao do estado da

arte sobre o tema com enfase nas aplicacoes encontradas na literatura. Em seguida,

apresenta-se de forma resumida os fundamentos basicos para a definicao formal de

acomodacao elastica e tambem os teoremas classicos que descrevem o fenomeno na

forma contınua e discreta. Posteriormente e fornecida em detalhes a formulacao de

tres elementos finitos: um elemento de viga bidimensional, um elemento de casca

axissimetrica e um elemento de casca tridimensional. A eficiencia dos elementos

propostos em analise de shakedown, e no caso particular de analise limite, e testada

atraves de um conjunto de exemplos numericos.

vii

Abstract of Thesis presented to COPPE/UFRJ as a partial fulfillment of the

requirements for the degree of Doctor of Science (D.Sc.)

CONTINUUM-BASED FINITE ELEMENTS FOR SHAKEDOWN ANALYSIS

OF PIPELINES AND PRESSURE VESSELS


September/2013

Advisors: Nestor Alberto Zouain Pereira

Lavinia Maria Sanabio Alves Borges

Department: Mechanical Engineering

In this work we present structural continuum-based (CB) elements devised for

the computation of collapse loads and shakedown limits of pipelines and pressure

vessels. Firstly, a review of the state of the art on the theme is provided especially

focusing on applications found in literature. Then, the basic knowledge for the

formal definition of shakedown is introduced as well as the classical shakedown

theorems in continuum and discrete setting. Next, formulations of three different

CB finite elements are presented in detail, namely: a two-dimensional beam element,

an axisymmetric shell element and, finally, a three-dimensional shell element. The

efficiency of the proposed elements in limit and shakedown problems is assessed by

means of selected numerical examples.

viii

Contents

List of Figures xii

List of Tables xv

1 Introduction 1

1.1 Research Objective . . . . . . . . . . . . . . . . . . . . . . . . . . . . 3

1.2 Literature review . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 4

1.3 Thesis Outline . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7

2 The Classical Shakedown Theory 9

2.1 Basic Notation and Definitions . . . . . . . . . . . . . . . . . . . . . . 10

2.1.1 Kinematics and equilibrium . . . . . . . . . . . . . . . . . . . 10

2.1.2 Elastic-ideally plastic materials . . . . . . . . . . . . . . . . . 10

2.2 Elastic Shakedown . . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2.1 Definition of elastic shakedown . . . . . . . . . . . . . . . . . 13

2.2.2 Load domain . . . . . . . . . . . . . . . . . . . . . . . . . . . 13

2.2.3 Conditions for shakedown . . . . . . . . . . . . . . . . . . . . 14

2.3 Extremum Principles for Elastic Shakedown . . . . . . . . . . . . . . 15

2.3.1 Statical maximum principle . . . . . . . . . . . . . . . . . . . 15

2.3.2 Mixed principles . . . . . . . . . . . . . . . . . . . . . . . . . 15

2.3.3 Kinematical minimum principle . . . . . . . . . . . . . . . . . 16

2.3.4 Optimality conditions in the continuum setting . . . . . . . . 16

2.4 Discrete formulation for elastic shakedown . . . . . . . . . . . . . . . 16

2.4.1 Discrete load domain . . . . . . . . . . . . . . . . . . . . . . . 17

2.4.2 Verification of plastic admissibility . . . . . . . . . . . . . . . 17

2.4.3 Compatibility and equilibrium . . . . . . . . . . . . . . . . . . 17

2.4.4 Extremum principles in the discrete setting . . . . . . . . . . . 18

3 Mixed Finite Elements in Linear Elasticity 21

3.1 Two field mixed formulation in elasticity . . . . . . . . . . . . . . . . 23

3.2 Discrete formulation in elasticity . . . . . . . . . . . . . . . . . . . . . 23

3.2.1 Interpolation functions . . . . . . . . . . . . . . . . . . . . . . 23

ix

3.2.2 Discrete equilibrium, compatibility and constitutive equations 24

3.2.3 Remarks on stability and convergence in elasticity . . . . . . . 25

4 Solution Algorithm for Shakedown Analysis 27

4.1 The linearization process . . . . . . . . . . . . . . . . . . . . . . . . . 27

4.2 The solution algorithm . . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.2.1 Newton iterations . . . . . . . . . . . . . . . . . . . . . . . . . 29

4.2.2 Relaxation and uniform scaling of the variables . . . . . . . . 31

5 Continuum-based Beam Element 33

5.1 Geometry, displacements and strains . . . . . . . . . . . . . . . . . . 33

5.1.1 Master and slave nodes . . . . . . . . . . . . . . . . . . . . . . 34

5.1.2 Kinematics of the underlying continuum element . . . . . . . . 36

5.1.3 The interpolation of displacements and velocities . . . . . . . 38

5.1.4 Enforcing bending theory hypotheses . . . . . . . . . . . . . . 39

5.1.5 Computing strain in a generic point . . . . . . . . . . . . . . . 40

5.2 Stresses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 41

5.3 Yield function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

5.4 Discrete strain operator . . . . . . . . . . . . . . . . . . . . . . . . . 43

5.5 Discrete elastic relation . . . . . . . . . . . . . . . . . . . . . . . . . . 46

5.6 Numerical Examples . . . . . . . . . . . . . . . . . . . . . . . . . . . 47

5.6.1 Linear elastic solution of a curved beam . . . . . . . . . . . . 48

5.6.2 Collapse of a beam in tension plus bending . . . . . . . . . . . 49

5.6.3 Shakedown analysis in tension plus bending . . . . . . . . . . 53

6 Continuum-based Axisymmetric Shell Element 56






6.1.5 Computing strain in a generic point . . . . . . . . . . . . . . . 61

6.2 Stresses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 62

6.3 Yield function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 63




6.6.1 Limit analysis of cylindrical shells under a ring load . . . . . . 66

6.6.2 Limit and shakedown analysis of pressure vessels with ellip-

soidal and torispherical heads . . . . . . . . . . . . . . . . . . 72

x

7 Continuum-based Three-Dimensional Shell Element 83






7.1.5 Computing strain at a generic point . . . . . . . . . . . . . . . 89

7.2 Stresses . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 91

7.3 Yield function . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 92



7.6 Eliminating singularities in rotational stiffness . . . . . . . . . . . . . 95


7.7.1 Straight pipe under combined loading . . . . . . . . . . . . . . 96

7.7.2 Shakedown analysis of pipe bend subjected to internal pres-

sure and bending moment . . . . . . . . . . . . . . . . . . . . 102

7.7.3 Limit analysis of a cylinder-cylinder intersection . . . . . . . . 105

7.7.4 Collapse of Dented Pipelines Subjected to External Pressure . 109

8 Conclusion 124

Bibliography 127

xi

List of Figures

1.1 Stress vs. strain response of structures subjected to cyclic load pro-

grams. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 2

1.2 Elastic-perfectly plastic material. Uniaxial test and yield surface in

principal stresses plane. . . . . . . . . . . . . . . . . . . . . . . . . . . 2

2.1 Associated plastic flow rule. (a) Differentiable yield surface. (b) Sub-

differentiable yield surface. . . . . . . . . . . . . . . . . . . . . . . . . 12

5.1 Continuum-based beam element . . . . . . . . . . . . . . . . . . . . . 34

5.2 Semi-circular arch subjected to a load 2P . . . . . . . . . . . . . . . . 48

5.3 Nodal displacements and rotation for the semi-circular arch conside-

ring various slenderness ratios . . . . . . . . . . . . . . . . . . . . . . 50

5.4 Deformation of the circular arch: analytical solution and FE results . 51

5.5 Beam under axial traction and bending moment . . . . . . . . . . . . 51

5.6 Collapse limits for a beam section under axial traction and bending

moment: analytical solution and FE results. Left: discontinuous

stress field between layers. Right: continuous stress field between

layers. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 52

5.7 Bree diagram for a beam under constant axial traction and a variable

uniform bending moment, ZOUAIN [10]. . . . . . . . . . . . . . . . . 54

5.8 Bree diagram for shakedown analysis of a beam under constant axial

traction and a variable uniform bending moment:analytical solution

and FE results. Left: discontinuous stress field between layers. Right:

continuous stress field between layers. . . . . . . . . . . . . . . . . . . 55

6.1 Continuum-based axisymmetric shell element: displacement nodes

(master and slave) and stress nodes . . . . . . . . . . . . . . . . . . . 57

6.2 Cylindrical shell under ring load . . . . . . . . . . . . . . . . . . . . . 66

6.3 Finite element meshes . . . . . . . . . . . . . . . . . . . . . . . . . . 69

6.4 Example of mesh refinement results: long cylinder (ξL = 3.536) . . . . 70

6.5 Collapse mechanisms of cylinders under ring load . . . . . . . . . . . 70

xii

6.6 Plastic strain rates for a short cylinder (ξL = 0.511) under collapse

due to a ring load . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

6.7 Collapse load versus non-dimensional cylinder’s length . . . . . . . . 72

6.8 Geometrical parameters of pressure vessels . . . . . . . . . . . . . . . 73

6.9 Collapse load versus thickness ratio for pressure vessels with ellipsoi-

dal heads . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 77

6.10 Examples of collapse mechanisms for pressure vessels with ellipsoidal

heads under internal pressure. . . . . . . . . . . . . . . . . . . . . . . 78

6.11 Examples of collapse mechanisms for pressure vessels with torisphe-

rical heads under internal pressure. . . . . . . . . . . . . . . . . . . . 79

6.12 Equivalent elastic stress plots from solid (left) and shell (right) models

for the pressure vessel with torispherical head. Geometrical parame-

ters of the vessel: a/b = 2.5, r/a = 0.12 and a/h = 10. . . . . . . . . 82

7.1 Continuum-based shell element: displacement nodes (master and

slave) and stress nodes . . . . . . . . . . . . . . . . . . . . . . . . . . 84

7.2 Finite element mesh for thin-walled pipe (r/t = 25) . . . . . . . . . . 97

7.3 Collapse load curve of thin-walled pipe (r/t = 25) under internal

pressure and bending . . . . . . . . . . . . . . . . . . . . . . . . . . . 97

7.4 Stress field in the cross section of a thick-walled pipe under maximum

combination of axial tension and bending . . . . . . . . . . . . . . . . 98

7.5 Collapse load curve of thick-walled pipe (r/t = 5) under axial force

and bending . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 99

7.6 Plastic strain rates for a thick-walled pipe (r/t = 5) under collapse

due to the loads (F/FC ,M/MC) = (0.4924, 0.7162) . . . . . . . . . . 100

7.7 Load domains for (a) the limit analysis and (b) the shakedown analy-

sis of a thick-walled pipe (r/t = 5) under internal pressure and bending100

7.8 Collapse load curve and shakedown limit for a thick-walled pipe

(r/t = 5) under internal pressure and bending . . . . . . . . . . . . . 101

7.9 Pipe bend geometry and loads . . . . . . . . . . . . . . . . . . . . . . 102

7.10 Finite element mesh using CB shell elements . . . . . . . . . . . . . . 103

7.11 Collapse load curve and shakedown limit for a thick-walled (r/t = 5)

and highly curved (h = 0.4) pipe bend under internal pressure and

bending. CL = collapse load; SD = shakedown limit . . . . . . . . . 104

7.12 Shell-solid model for the pipe bend. Detail of the coarser mesh em-

ployed in bend region . . . . . . . . . . . . . . . . . . . . . . . . . . . 105

7.13 Equivalent elastic stress (von Mises) from ANSYS mesh 3 . . . . . . . 106

7.14 Equivalent elastic stress (von Mises) from CB shell model . . . . . . . 106

7.15 Cylinder-cylinder intersection: geometry and loads . . . . . . . . . . . 107

xiii

7.16 Cylinder-cylinder intersection: finite element mesh with CB shell ele-

ments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

7.17 Cylinder-cylinder intersection: finite element mesh with ANSYS

SHELL93 elements . . . . . . . . . . . . . . . . . . . . . . . . . . . . 108

7.18 Collapse load curve for the cylinder-cylinder intersection . . . . . . . 109

7.19 Geometric parameters of most deformed cross section of a dent . . . . 113

7.20 Geometry of the dented pipe model . . . . . . . . . . . . . . . . . . . 116

7.21 Boundary Conditions applied in 3-D models . . . . . . . . . . . . . . 117

7.22 Coarser mesh used in the dented pipe models . . . . . . . . . . . . . 118

7.23 Finer mesh used in the dented pipe models . . . . . . . . . . . . . . . 119

7.24 Pressure vs. displacement curve for pipes with R/t = 9.45. . . . . . . 119

7.25 Example of plastic strain rate distribution for the pipe with a dent

depth dSB/D = 24. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 121

7.26 Maximum pressure vs. ovalization. . . . . . . . . . . . . . . . . . . . 122

7.27 Deformed shape of a pipe with an initial ovalization of 0.5% for dif-

ferent levels of external pressure. . . . . . . . . . . . . . . . . . . . . 123

xiv

List of Tables

5.1 Maximum differences between analytical and FE results . . . . . . . . 53

6.1 Load multipliers (µ) obtained with CB shell models and solid axisym-

metric models. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71

6.2 Limit loads calculated with CB shell element for pressure vessels with

ellipsoidal heads. pE = elastic limit; pLA = collapse load; pSD =

shakedown load; Mch = failure mechanism; AP = alternate plasticity;

IC = incremental collapse. . . . . . . . . . . . . . . . . . . . . . . . . 75

6.3 Limit loads calculated with CB shell element for pressure vessels with

torispherical heads. pE = elastic limit; pLA = collapse load; pSD =

shakedown load; Mch = failure mechanism; AP = alternate plasticity;

IC = incremental collapse. . . . . . . . . . . . . . . . . . . . . . . . . 76

6.4 Difference between limit loads calculated with CB shell and results

from YEOM and ROBINSON [76] for pressure vessels with ellipsoidal

heads. pE = elastic limit; pLA = collapse load; pSD = shakedown load 79

6.5 Difference between limit loads calculated with CB shell and results

from YEOM and ROBINSON [76] and MAKRODIMOPOULOS [72]

for pressure vessels with torispherical heads. pE = elastic limit; pLA

= collapse load; pSD = shakedown load . . . . . . . . . . . . . . . . . 80

7.1 Load multiplier µ preventing alternate plasticity for the pipe bend

under cyclic bending moment . . . . . . . . . . . . . . . . . . . . . . 107

7.2 Non-dimensional geometric parameters of the dented pipe model . . . 117

7.3 Collapse pressures calculated by incremental analysis. . . . . . . . . . 120

7.4 Collapse pressures calculated by direct analysis. . . . . . . . . . . . . 121

7.5 Collapse pressures calculated by direct analysis considering the final

geometry from the incremental analysis . . . . . . . . . . . . . . . . . 121

xv

Chapter 1

Introduction

The safety assessment of structures subjected to variations of loading is a common

and very important task of structural analysis. Besides the verification of possible

failure due to high cycle fatigue, international codes, e.g. [1–3], require that load

fluctuations must not produce in the structure excessive accumulated plastic strains

which, in turn, can lead to either low cycle fatigue failure (alternate plasticity) or

unserviceability due to progressive incremental inelastic deformation (ratcheting).

Verification of ratcheting is generally mandatory in the codes whereas the determi-

nation of bounds for loading preventing low cycle fatigue is optional.

Shakedown of a structure occurs if, after a few cycles of load application, rat-

cheting ceases. The subsequent structural response is elastic, or elastic-plastic, and

progressive incremental inelastic deformation is absent. Elastic shakedown is the

case in which the consecutive response is elastic, which means elastic adaptation.

The different types of stress-strain behavior presented by a structure under cyclic

loads is illustrated in Fig 1.1.

Elastic shakedown of a structure can commonly be assessed by means of sha-

kedown analysis. In the classical shakedown theory [4–7] the structure material is

assumed linear elastic–perfectly plastic (ideal plasticity). For this type of material

model there is no hardening and the yield stress level is independent on plastification.

In other words, the yield surface remains fixed for any deformation process suffe-

red by the material. A example of a typical uniaxial cyclic (tension–compression)

stress–strain curve for a elastic-perfectly plastic material is shown in Fig. 1.2(a).

Figure 1.2(b) depicts the corresponding representation of the yield surface for a von

Mises model in the space of principal stresses [8].

For structural components made of ideally plastic material (or any material pre-

senting limited work-hardening), it is theoretically possible to determine the limit

load that produces instantaneous collapse of the structure. In this case, instan-

taneous collapse takes place when a non-fluctuating load produces kinematically

admissible plastic strain rates under constant stresses. This phenomenon can be

1

Figure 1.1: Stress vs. strain response of structures subjected to cyclic load programs.

Figure 1.2: Elastic-perfectly plastic material. Uniaxial test and yield surface inprincipal stresses plane.

2

seen as a particular case of ratcheting (for zero amplitude of load fluctuations) and

can be also treated within the classical shakedown framework [9, 10]. In this work

we employ the shakedown theory to calculate both collapse and shakedown limit

loads.

We remark that the classical shakedown analysis has some important advanta-

ges as compared to classical incremental plasticity procedures: the safety factor that

guarantees elastic adaptation is directly computed without resorting to theoretically

infinite number of full incremental load-displacement analyses. Due to this reason

the classical shakedown analysis is called a direct method. Moreover, in the classi-

cal shakedown theory the safety assessment of the structure is made with the sole

knowledge of the load range rather than full load history. Therefore, the shakedown

analysis allows the engineer to work with more realistic hypotheses in terms of loads.

The assessment procedure in the classical shakedown theory generally involves

finite element discretizations of the structure and the use of dedicated numerical

procedures for the solution of the shakedown analysis problem [9]. In the case of

thin structures, such as pipes and pressure vessels, three-dimensional solid finite

element models can be computationally very expensive due to the large number

of elements required. Coarser meshes, on the other hand, are prone to contain

elements with high aspect ratios which in turn can cause ill conditioning of the

finite element equation system and lead to loss of accuracy in the solution. These

problems are generally circumvented by using structural elements such as beam,

plate and shell elements. However, for limit and shakedown analysis, the use of

classical structural elements requires the development of a solution approach with

laborious approximations in terms of generalized stresses to express, for instance,

the yield limit of the material [9].

1.1 Research Objective

The objective of this work is the development of structural elements for collapse

and shakedown calculation which, in particular, bypass the need to express the

shakedown formulation in terms of generalized stresses. Then, the problem is posed

in a simpler manner and one can easily combine the elements with any shakedown

algorithm developed for classic solid elements without major modifications.

The proposed elements are formulated using the continuum-based (CB) approach

[11–18] and a mixed Hellinger-Reissner-type formulation [19–21] with stress and dis-

placement/velocity fields interpolated independently and, to the best our knowledge,

this is the first report proposing this type of structural elements for applications in

collapse and shakedown analysis.

In the author’s opinion the use of mixed elements in this context is advantageous

3

because they can be more easily implemented. In addition, for kinematically-based

formulations [12] the use of the kinematic minimum principle for shakedown analysis

requires the computation of the dissipation function, which is unbounded unless

plastic strain rates are plastically admissible at all point of the element [21]. This

introduces a very stringent constraint in non-constant plastic strain interpolations.

In the present development the attention is focused on the ability of the finite

elements to efficiently represent solutions of the considered problem, with special

emphasis on the accurate determination of critical amplification factors for the pres-

cribed reference domain of variable loadings. The numerical procedure used here for

the general shakedown problem is presented in detail in [10, 22].

1.2 Literature review

According to MAIER et al. [9] the history of direct methods can be dated back to

the mid-17th century. They cite Galileo’s last “dialog” as the first work where a

kinematic approach for limit analysis is formulated [23].

In 1932 Bleich states the first elastic adaptation theorem for hyperstatic beams,

made of elastic-perfectly plastic material [6, 10, 24]. Bleich’s theorem, in turn, has

its origin in an work of Gruning, published three years before, where the static

shakedown theorem was proved for a system of beams with ideal I cross section

[25, 26].

In 1938 Melan generalizes Bleich’s theorem for the three-dimensional continuum

and provides theorems for two different types of material: elastic-perfectly plastic

solids and solids presenting limited linear kinematic hardening [27, 28]. In 1956

Koiter states the kinematic theorem for shakedown under the hypothesis of elastic-

perfectly plastic material [29]. The dual theorems of Melan’s and Koiter’s are the

basis of the classical shakedown theory.

Besides the assumption of ideal plasticity, Melan’s theorem is based on the fol-

lowing simplifying hypotheses [25]:

• Infinitesimal displacements and strains;

• Isothermal and quasi-static process;

• Time independent plasticity with associated plastic flow rule and a convex

yield surface in stress space;

• The free-stress state is strictly interior to the elastic domain;

• Mechanical properties independent from temperature;

• All parameters are deterministic quantities;

4

• The only goal of the analysis is the computation of a safety factor making sure

elastic adaptation.

Several works have been developed from the pioneer theorems of Melan and

Koiter with the aim of providing formulations under less stringent assumptions,

and thus expanding the range of application of the shakedown theory. As a result,

there are in literature extensions of the classical shakedown theory that can account

(in some degree) for work-hardening, temperature effects, non-associated plasticity,

inertia and damping, finite strains and finite displacements, creeping, damage etc.

The works of KONIG and MAIER [30], WEICHERT and MAIER [31], MROZ et al.

[32], PYCKO and MAIER [33], NGUYEN [34], ZOUAIN [10], MAIER et al. [9] and

ABDEL-KARIM [26] provide a good picture of the evolution of shakedown theory

and list some important contributions.

It is worthy to mention the efforts that have been made in the last two decades

to implement algorithms of direct methods for shakedown and limit analysis in

commercial finite element programs such as ANSYS, ABAQUS, ADINA, and so

forth [25]. With such implementations the capabilities for pre- and post-processing

as well as the vast libraries of finite elements existing in commercial programs could

be used for larger applications of shakedown analysis with industrial purposes.

Recently, in the period from 1998 to 2002, it was carried out in Europe the project

“FEM-Based Limit and Shakedown Analysis for Design and Integrity Assessment

in European Industry – LISA”. The objective of that project was to combine direct

methods of plasticity with stochastic analysis by means of First and Second Order

Reliability Methods (FORM/SORM) in a general purpose Finite Element package.

The results of LISA project are available in the report edited by STAAT et al. [35]

and in scientific publications, such as [36–41].

Limit and shakedown analysis have a significant number of applications in struc-

tural analysis. We classify the numerical examples found in literature in two ca-

tegories. The first class of examples comprises structures of simpler geometry and

loads, usually two-dimensional or axisymmetric. Most of these problems has known

analytical solution. They are especially devised as benchmarks to approximate te-

chniques, which generally make use of the finite element method. The majority of

published numerical examples found in papers falls into this category. Illustratively

we refer to the following examples in shakedown analysis:

• Clamped beams under constant axial load plus variable load (bending moment

or transverse load) applied at its free end [42, 43];

• Perforated plates subjected to variable tension [22, 25, 42–46];

• Thick-walled pipes under variable internal pressure [42, 47];

5

• The Bree-problem: thin-walled pipe subjected to fixed internal pressure and

variable linear temperature gradient across its wall [46, 48];

• Pressure vessels with torispherical head under variable internal pressure [25,

36].

With the development of large-scale optimization algorithms, a second class of

examples comprising more complex geometries and loads began to appear in research

papers, such as:

• Shakedown analysis of pipe bends under variable internal pressure, bending

moment and/or temperature [49–51];

• Shakedown analysis of pipes containing single or multiple defects (wall loss)

subjected to variable axial force, bending moment and internal pressure [38,

52–54];

• Limit analysis of pipes with cracks under internal pressure [55, 56];

• Shakedown analysis of pipe junctions under variable internal pressure and

temperature [36, 39];

Due to their relevance, some examples above, although not necessarily new, have

been revisited and solved with up-to-date shakedown algorithms and/or elements

presenting a better performance. The results found in the more recent papers are

usually compared with available simplified analytical solutions and previous nume-

rical results.

Three-dimensional models found in literature are very frequently built with solid

finite elements, although some geometries could be more efficiently discretized in

terms of computational effort by using shell elements. However, it is not easy to

find in research papers and reports examples where shell elements are employed in

limit analysis and especially in shakedown analysis. A remarkable exception are the

axisymmetric structures which have been solved using shell formulations for a long

time [57–59]. A possible reason to avoiding the use of shell elements in shakedown

formulations is the complexity inherent to the definition of admissibility of stresses

in terms of shell forces. We readily overcome this problem with the CB elements

proposed in this thesis. This is one of the motivations for the present work.

In addition, the good accuracy exhibited by the family of mixed stress-velocity

finite elements employed by Zouain and collaborators in limit and shakedown analy-

sis [10, 60] stimulated the development of structural finite elements with the same

interpolation scheme.

Both the continuum-based approach [11–18] and the two-field mixed interpola-

tion technique [19, 61] are not novel in the finite element method; however we could

6

not find previous works in literature where structural finite elements similar to that

presently proposed were applied in shakedown analysis or even in limit analysis. In-

deed, it is difficult to find works on the development of CB shell elements with any

mixed interpolation scheme for limit and shakedown analysis. For three-dimensional

applications, such elements are much rarer. In this review we could only find the

works of Corradi and Panzeri [62, 63] who devised three-dimensional shell elements

employing a different interpolation scheme for application in sequential limit analy-

sis.

The combination of the CB approach with the aforementioned mixed stress-

velocity interpolation strategy showed simple and effective for the examples discus-

sed in the following chapters.

1.3 Thesis Outline

This work is organized in two parts. The theoretical framework employed in the

thesis is presented in the first part which comprises chapters 2 to 4.

We begin in Chapter 2 with the presentation of the basic notation adopted to

address ideal plasticity, then we summarize the classical shakedown theorems which

are given in the form of extremum principles in the continuum and discrete setting.

Next, in Chapter 3 the formulation of mixed stress-displacements finite elements

in the context of elasticity is briefly discussed. The objective is to present the

regularity required for the trial functions according to the respective variational

principle, i.e., the Hellinger-Reissner functional. Moreover, we outline some stability

and convergence issues for this type of finite element.

Then, in Chapter 4 the algorithm employed for the solution of the shakedown

problem is presented in brief. Details are omitted because the elements presented

are applicable to other solution strategies and, more importantly, because it is our

intention to focus on the formulation of the finite elements and in their applications,

discussed in the chapters that follows.

The second and last part of the thesis comprises the chapters 5 to 8. This

part comprises the detailed description of the finite elements proposed in this work

and examples of application that confirm the accuracy of the finite elements in the

computation of collapse loads and shakedown limits.

This research was planned and executed so that the elements could be formulated

and tested in stages with increasing level of complexity. In this text the order of

the chapters reflects the way that the work was carried out. Thus, in Chapter 5 the

CB beam element is presented and tested in a small, but sufficient, set of examples

that gave us confidence to go to the next step: the development of the axisymmetric

shell element.

7

The CB axisymmetric shell element presented in Chapter 6 is easily derived from

the CB beam element. The CB axisymmetric shell element is more extensively tested

in limit and shakedown analysis with two representative examples comprising pipes

and pressure vessels. The results obtained with this proposed element are compared

with numerical and analytical results from literature and also with outcomes from

incremental analyses performed in a commercial finite element code.

In a similar way, in Chapter 7 the CB three-dimensional shell element is pre-

sented. This element is tested in collapse and shakedown analysis of straight and

curved pipes as well as in the computation of collapse loads for a pipe junction and

dented pipes. Again results obtained with the proposed CB shell element are verified

against both analytical and numerical results.

Finally the work is closed with the conclusions listed in Chapter 8.

8

Chapter 2

The Classical Shakedown Theory

In this chapter we provide a brief review of the classical shakedown theory [4, 6],

mainly following [10]. In particular we present here: (i) the formal definition of sha-

kedown; (ii) the basic theorems of the classical theory; and (iii) the set of equations

to be solved in shakedown problems both in the continuum and the discrete setting.

Extensions of shakedown theory for taking into account hardening and temperature

effects, dynamic effects, geometrical nonlinearities, etc. and for the identification of

the collapse mechanism, also addressed in [10], are not discussed in this work.

After the introduction of some important basic notation and definitions, the

elastic shakedown is formally given. Then, classical theorems of Melan and Koiter,

which determine the conditions for shakedown to take place in a continuum body

subjected to variable loads, are presented. Further, the theorems aforementioned

are rewritten as extremum principles, named: statical, mixed and kinematic. Fi-

nally, the continuum formulation is discretized and provided as a set of optimality

conditions which can be solved using numerical procedures.

The classical shakedown theory deals with the prevention of three possible modes

of failure presented by structures made of perfectly plastic material when subjec-

ted to the combination of fixed and variable loading, namely: alternating plasticity

(plastic shakedown or low cycle fatigue), incremental collapse (ratcheting), or ins-

tantaneous collapse (plastic collapse).

The main objective of shakedown analysis is the computation of the load amplifi-

cation factor, µ , ensuring elastic adaptation. The computation of the amplification

factor that only prevents instantaneous collapse is particularly called limit analysis.

Shakedown theory allows working under the realistic assumption that only the

range of variable loading is known. The procedures developed under this theory,

collectively named direct methods, make it possible the direct computation of the

amplification factor without the need of a theoretically generally infinite number of

full incremental analyses. Moreover, only the theory of direct methods can answer

whether critical loads or cycles do exist, or not, independently from load histories.

9

2.1 Basic Notation and Definitions

In this section the notation used for the exposition of the classical theory of sha-

kedown is briefly provided. Some basic concepts of convex analysis [64] that are

usually applied in plasticity are also presented.

2.1.1 Kinematics and equilibrium

Consider a body occupying an open bounded region B with regular boundary Γ.

The set of all admissible velocity fields v complying with homogeneous boundary

conditions prescribed on the part Γu of Γ is the space V . Strain rate tensor fields

d are elements of W , and the tangent deformation operator D maps V into W .

The dual space of W is the space W ′ of stress fields σ. The equilibrium operator

D′, dual of D, maps elements of W ′ into the space of load systems F denoted by

V ′. Prescribed loads F vary quasi-statically. The spaces V and V ′ are also dual.

Consequently, static and kinematic relations are written as

d = Dv, F = D′σ. (2.1)

Moreover, the hypothesis of small deformation holds and the infinitesimal displace-

ment and strain fields, u and ε respectively, are related by

ε = Du. (2.2)

The internal power for any pair σ ∈ W ′ and d ∈ W is given by the duality

product ⟨σ, d

⟩:=

∫Bσ(x) · d(x)dB (2.3)

and the equation of external power for any pair F ∈ V ′ and v ∈ V is⟨F , v

⟩:=

∫Bb · vdB +

∫Γτ

τ · vdΓτ , (2.4)

where b and τ represent volume and boundary load densities.

We also define the set of residual (self-equilibrated) stress field Sr, such that

Sr :=σ |⟨σ, d

⟩= 0∀d ∈ W (2.5)

2.1.2 Elastic-ideally plastic materials

The additive decomposition of total strain ε into an elastic εe and an inelastic εp

part is assumed. The general specific and global Helmoltz free energy are denoted

10

by W(ε(x)) and W(ε) =∫BW(ε(x))dB respectively. Particularly, the following

expressions can be written when the linear elastic relation is assumed

ε = εe + εp, σ = Eεe (2.6)

and

W(εe) =1

2

⟨Eεe, εe

⟩, Wc(σ) =

1

2

⟨σ, E−1σ

⟩, (2.7)

where E denotes the elastic operator and Wc the complementary free energy.

In addition, for any point x of an elastic-plastic body B, the stress tensor σ is

such that the plastic admissibility condition is fulfilled. In general, this condition is

imposed by means of m constraints, or yield modes; however, in this work we only

consider the von Mises yield criterion (i.e. m = 1). Therefore σ has to belong to

the set

P = σ | f(σ(x)) ≤ 0 ∀x ∈ B, (2.8)

where f is a real valued regular convex function of σ describing the yield criterion.

Additionally, we assume that the stress-free state is strictly admissible, i.e. σ(x) =

0 ∈ P .

Likewise, the closed convex set P of plastically admissible stress fields is defined

by

σ ∈ P⇔ σ(x) ∈ P ∀x ∈ B. (2.9)

By defining the plastic dissipation function per unit of volume as

D(dp(x)) = supσ∗(x)∈P

(σ∗(x) · dp(x)) (2.10)

and the indicator function IP (σ(x)) of P , that equals zero for any σ(x) ∈ P and

+∞ otherwise, then the constitutive relation between plastic strain rates dp(x) and

stresses σ(x) is written, for the case of associative plastic flow, as

σ(x) ∈ ∂D(dp(x))⇔ dp(x) ∈ NP (σ(x)) (2.11)

where the subdifferential ∂D(dp(x)) is the set of all stress tensors σ(x) such that

D(dp∗)−D(dp(x)) ≥ σ(x) · (dp

∗− dp(x)) ∀dp

∗(2.12)

and NP (σ(x)) := ∂IP (σ(x)) is the cone of normals to P at σ(x), i.e. the set of all

plastic strain rates dp(x) such that

(σ(x)− σ∗) · dp(x) ≥ 0 ∀σ∗ ∈ P. (2.13)

11

Figure 2.1: Associated plastic flow rule. (a) Differentiable yield surface. (b) Sub-differentiable yield surface.

The dissipation function D(dp(x)) can be identified as the support function of

P , hence it is sublinear, i.e. convex and positively homogeneous of degree one. It

is also non-negative, since σ(x) = 0 ∈ P . Therefore, it has minimum equal to zero,

attained at dp(x) = 0.

The model described above belongs to the class of standard materials presented

by HALPH and NGUYEN [65] apud [10]. Thus the second law of thermodynamics

is automatically fulfilled, i.e.

dint = σ(x) · dp(x) ≥ 0 (2.14)

with dint denoting the specific internal dissipation.

The material relations (2.11) are equivalent to the following classical form valid

for unimodal yield functions like the von Mises criterion: the plastic strain rate is

related to the stress, at any point of B, by the normality rule

dp(x) = λ(x)∇f(σ(x)). (2.15)

Here ∇f denotes the gradient of f with respect to σ, and λ is the field of plas-

tic multipliers. At any point of B, the plastic factors λ are related to f by the

complementary condition

λ ≥ 0, f ≤ 0 and f(σ(x)) · λ(x) = 0. (2.16)

Figure 2.1 graphically illustrates the plastic strain rates for a smooth yield func-

tion such as the von Mises model and for a more general and non-smooth yield

surface, where dp ∈ NP .

Note that, for the sake of simplicity of notation, the explicit dependence on

the point x ∈ B will be dropped whenever there is no possible misunderstanding

12

concerning a field, such as σ and its value σ(x).

2.2 Elastic Shakedown

The objective of this section is to formally define shakedown and present the neces-

sary and sufficient conditions for this phenomenon to happen.

2.2.1 Definition of elastic shakedown

Elastic shakedown takes place when the total plastic work, Wp(t) =∫ t

0

∫B dintdBdt,

is limited, i.e.

limt→∞

Wp(t) <∞ (2.17)

and, additionally, when changes in the residual stress field vanishes after a arbitrarily

large period in time during the load cycles. Hence,

∃ limt→∞

[σ(t)− σE(t)] (2.18)

In the equation above σ(t) denotes the stress field in the actual body produced by

a variable loading. The symbol σE(t) represents the elastic (unlimited) stress field

produced by the same loading in the reference body, which is identical to the real

one with respect to elastic behavior but is free of plastic admissibility constraints

on the stresses.

2.2.2 Load domain

An important feature in the formulation of shakedown theory is that loads are not

necessarily described as functions of time but rather assumed to arbitrarily vary

in a prescribed domain. This load domain is then mapped, by a purely elastic

transformation, to a set ∆E ⊂ W ′ of stress fields.

In what follows we define the pointwise envelope ∆ of the domain of elastic

stresses ∆E. This set ∆ is the preferred form for loading data in shakedown theory.

The set of local values of elastic stresses associated to any feasible loading is

∀x ∈ B ∆(x) := σE(x) | σE ∈ ∆E, (2.19)

then, the pointwise envelope of the set ∆E is defined as

∆ :=σ | σ(x) ∈ ∆(x) ∀x ∈ B

⊃ ∆E. (2.20)

13

Likewise, a superposed fixed (non-amplifying) load is represented by the corres-

ponding elastic stress field σF .

2.2.3 Conditions for shakedown

In this subsection the formulation of conditions for shakedown assessment in terms

of stresses are presented. This formulation, often known as Melan´s theorem, is

based on the works of BLEICH in 1932 [6] and MELAN [27].

Sufficient condition for shakedown

Melan-Koiter´s theorem provides the sufficient condition for shakedown to take place

in the sense that the total plastic work Wp is bounded irrespective to initial condi-

tions. It states that: “an elastic ideally plastic body will shakedown under a given

loading history if there exists a time-independent residual stress field σr, a scalar

m > 1 and an instant t0 such that the fictitious elastic stresses σE, produced by the

loading program acting on the unlimitedly elastic reference body, when amplified by

the factor m and superposed to the fixed residual stress field σr render plastically

admissible stresses ever after t0 ”. That is, if there exists σr ∈ Sr, m > 1 and t0

such that for all t ≥ t0

mσE(t) + σr ∈ P (2.21)

the body accommodates with Wp <∞.

The theorem above can be also extended to loading histories pertaining to the

prescribed loading domain ∆. In this case, elastic shakedown will occur if there

exist σr ∈ Sr and m > 1 such that

m∆ + σr ⊂ P (2.22)

Necessary conditions for shakedown

Koiter´s theorem complements Melan´s theorem and gives the non-shakedown con-

ditions. It states that “shakedown is impossible if no time-independent distribution

of residual stresses can be found with the property that the sum of residual stresses

and “elastic” stresses is an allowable state of stress (i.e. non strictly interior to P )

at every point of the body and for all possible load combinations”.

For a load domain, the theorem above is rewritten as: If, for some m > 1, there

is no residual stress field σr ∈ Sr such that σr +m∆ ⊂ P, the body does not shake

down to the (amplified) set m∆ of possible loads.

14

Extension of the shakedown theorem to consider superposition of fixed

loading to the variable loading

The shakedown theory is extended for taking into account the superposition of

variable loading to fixed loads. The fixed loads are represented by the stress field

σF obtained from the elastic response of the reference body subjected to the fixed

loading.

Thus, elastic shakedown is ensured for the prescribed domain of load variations,

denoted by ∆, in the presence of a fixed (non-amplified) load, represented by this

associated elastic stress field σF , if there exists a fixed residual stress field σr ∈ Sr

and a scalar m > 1 such that

σF +m∆ + σr ⊂ P (2.23)

2.3 Extremum Principles for Elastic Shakedown

In this section the shakedown theory presented in the previous section is rewritten

as extremum principles. The objective is to compute the limit factor µ that ensures

elastic shakedown for any variable stress histories belonging to µ∆.

2.3.1 Statical maximum principle

The shakedown factor is obtained by Bleich and Melan’s statical formulation [6] as

µ = sup(µ∗,σr)∈R+×W ′

µ∗ | σF + µ∗∆ + σr ⊂ P; σr ∈ Sr

, (2.24)

where R+ is the set of nonnegative real numbers.

In geometrical terms, solving (2.24) is equivalent to finding the translation

σr ∈ Sr in the stress space that moves the prescribed set ∆ allowing its maxi-

mum amplification. The resulting shifted and scaled domain is entirely contained

in the plastically admissible set P.

2.3.2 Mixed principles

The primal mixed principle is derived from the statical formulation by introducing

the self-equilibrium constraint as a penalty function, resulting in

µ = sup(µ∗,σ)∈R+×W ′

infv∈V

µ∗ + 〈σ,Dv〉 | σF + µ∗∆ + σ ⊂ P

. (2.25)

The dual mixed principle is achieved by formally interchanging the sup and

15

inf operations. This procedure, in general, does not guarantee that primal and dual

problems have the same solution. Indeed, some conditions must be verified to assure

this equivalence [64].

2.3.3 Kinematical minimum principle

The kinematical principle can be obtained from the dual mixed formulation [10]. In

this case,

µ = infv∈V

µ(Dv) (2.26)

where

µ(d) = sup(µ∗,σ)∈R+×W ′

µ∗ +⟨σ, d

⟩| σF + µ∗∆ + σ ⊂ P (2.27)

2.3.4 Optimality conditions in the continuum setting

We can write the same problem as the following set of optimality equations [10]:

σr ∈ Sr, (2.28)

〈σ,Dv〉 = 1, (2.29)

and for all x ∈ Bdp = Dv, (2.30)

dp ∈ NP (σF + µσ + σr), σ ∈ ∆. (2.31)

The above optimality conditions can be interpreted in terms of mechanics as

follows: Eq. (2.28) is the self equilibrium constraint; Eq. (2.29) imposes that the

standard solution dissipates a power equal to one; Eq. (2.30) implies that the strain

rates Dv are purely plastic; and Eq. (2.31) guarantees plastic admissibility of strains.

2.4 Discrete formulation for elastic shakedown

In this section the shakedown formulations are discretized for obtaining approximate

solutions by numerical methods. To this end we employ the finite element method.

The symbols used for fields in the continuum model are employed from now on

for their discrete counterparts without the superposed hat (ˆ).

The general steps to obtain the discrete formulations for shakedown are:

1. The load domain is given by a finite number n∆ of extreme loads. Thus, the

local domain of variable loading ∆(x), for all x ∈ B, is a convex polyhedron;

16

2. The continuum body B is divided in nel parts, producing a finite dimensional

model;

3. Plastic admissibility constraints are verified in a set xj; j=1:p of critical

points (vertices) judiciously chosen in the body according to the discretization

(the notation 1:p abbreviates the sequence of integers from 1 to p).

2.4.1 Discrete load domain

After discretization, the load domain is given by the convex hull of the set of extreme

loadings

∆ = coσk; k=1:m

, m := p n∆. (2.32)

2.4.2 Verification of plastic admissibility

The plastic admissibility condition σF + µ∆ + σr ⊂ P, verified for each m = pn∆

local vertex of the mesh, can be given in the equivalent form

σF + µσk + σr ∈ Pk, k=1:m, (2.33)

where Pk represents the elastic range of a particular point in the body, written in

terms of the global parameters, i.e. only the pertinent components are selected in

advance.

For the continuum-based finite elements described in this work plastic admissibi-

lity constraints are verified at the stress nodes, which are placed at each corner of the

element layers. The stress fields are interpolated in the element by means of stan-

dard Lagrange linear functions and, consequently, plastic admissibility is assured for

any point within the layers. Under the hypothesis of homogeneous material in each

element layer, the property of plastic admissibility follows from: (i) the convexity of

the plastically admissible stress domain, and (ii) due to the stress at internal points

in each element layer to be convex combination of the values of the layer vertices.

2.4.3 Compatibility and equilibrium

The discrete equations for compatibility and equilibrium are respectively

d = Bv, (2.34)

and

BTσr = 0, (2.35)

17

where the global strain rate and residual stress vectors are denoted respectively by

d ∈ Rq and σr ∈ Rq. The symbol B : Rn → Rq denotes the discrete strain operator.

Superscripts n and q respectively denote the number of degrees of freedom that

remains after excluding rigid body motions and the number of stress components

for the entire mesh. The symbol BT is the discrete equilibrium operator which is

the adjoint of B. Note that in the above, and from now on, we use the superscript

T to denote transpose of a matrix.

It is worth emphasizing that the final form of B and, consequently, BT depends

on the discretization approach and the variational principle used. The respective

vectors of residual stress and strain rate must be interpreted accordingly. For ins-

tance, in this work we employ an interpolation strategy where both velocity and

stress fields are approximated. Interpolated fields are given, in this case, by the

following general form:

v(x) ≈ v(x) = Nvvi (2.36)

σ(x) ≈ σ(x) = Nσσi (2.37)

where: Nv(x) and Nσ(x) denote the interpolation functions for velocity and stress

fields respectively, and vi and σi are the vectors which collect the displacement and

stress nodal coefficients.

In the discrete mixed formulation the discrete strain operator of a finite element

i is obtained by replacing (2.36) and (2.37) in (2.25). Thus,

Bi :=

∫BiNTσDNvdBi. (2.38)

The matrix B is assembled with standard procedures of the finite element method.

Then, considering form the discrete strain operator (2.38), obtained for the discrete

mixed formulation, and in view of (2.34), one can recognize the nodal strain rate

component vector di ∈ Rq as collecting approximations of the corresponding compa-

tible strains at each vertex multiplied by the element volume. Note, however, that

the form of the operator Bi, which comprises interpolation functions for stress and

velocities, would be different if pure statical or kinematical formulation is used. The

meaning for the vectors of parameters, such as the vector di, would be different as

well.

2.4.4 Extremum principles in the discrete setting

As a consequence of the discretization steps, the extremum principles provided in

section 2.3 can be now rewritten in the following discrete forms [10]:

18

• Discrete statical formulation:

µ = supµ∗∈R+

σr∈Rq

µ∗ | σF + µ∗σk + σr ∈ Pk, k=1:m; BTσr = 0

(2.39)

• Discrete mixed formulation:

µ = supµ∗∈R+

σ∈Rq

infv∈Rn

µ∗ + σ ·Bv | σF + µ∗σk + σ ∈ Pk, k=1:m

(2.40)

• Discrete kinematic formulation:

µ = infv∈Rn

dk∈(Rq)m

∑Dk(dk)−σF ·Bv |

∑dk = Bv;

∑σk ·dk = 1

(2.41)

where

Dk(dk) = supσ∗∈Pk

σ∗ · dk (2.42)

Finally, the set of relations characterizing the solution of the problems (2.39-

2.41) is considered. These relations are properly derived using the corresponding

Lagrangian function

L(v,σr, µ, dk) = µ+ σr ·Bv −∑

[(σF + µσk + σr) · dk −Dk(dk)] (2.43)

Then, invoking the stationarity of L relative to the variables v, µ and σr yields

∂L∂v

= 0⇒BTσr = 0 (2.44)

∂L∂µ

= 0⇒∑

σk · dk = 1 (2.45)

∂L∂σr

= 0⇒∑

dk = Bv (2.46)

Equations (2.44) to (2.46) are the respective discrete counterparts of the equilibrium

constraint (2.28), the standard dissipated power (2.29) and compatibility between

plastic strain rate and velocity fields (2.30).

The stationarity of L in respect to dk is not as straightforward as presented

above for the other arguments of this function. Note that the dissipation function

Dk is not differentiable, thus, sub-differentials of L with respect to dk are calculated

instead of partial derivatives. Then, by taking into account that 0 ∈ ∂D, theorems

of convex analysis can be applied resulting in the following expression

∂L∂d3 0⇒ σF + µσk + σr ∈ ∂Dk(dk) k=1:m. (2.47)

19

The equation (2.47) is the discrete counterpart of (2.31). It states that the flow

produced by the stresses σF +µσk+σr is the strain rate dk. Flow takes place during

the critical load cycle when the stress associated to the load vertex k is active (i.e.

in the yield surface). Consequently, the plastic strain rate dk may produce a non-

vanishing plastic flow to the admissible plastic strain rate cycle Bv, as stated in

(2.46).

For the numerical algorithm used in the solution of the shakedown problem it is

particularly convenient rewriting the above set of optimality conditions (2.44)-(2.47)

by taking into account that the von Mises yield criterion in the form f(σ) ≤ 0 is a

scalar regular convex function. Hence, strain rates and the yield function gradient

are related by dk = λk∇f(σF + µσk + σr). Thus, the optimality conditions can be

equivalently stated as

• Discrete optimality conditions for unimodal yield function:

BTσr = 0 (2.48)∑k=1:m

λk[σk · ∇f(σF + µσk + σr)] = 1 (2.49)∑

k=1:m

λk∇f(σF + µσk + σr) = Bv (2.50)

λkf(σF + µσk + σr) = 0, k=1:m (2.51)

f(σF + µσk + σr) ≤ 0, k=1:m (2.52)

λk ≥ 0, k=1:m. (2.53)

20

Chapter 3

Mixed Finite Elements in Linear

Elasticity

In linear elasticity, mixed finite elements can be obtained from variational principles

that can be viewed as extensions of the principle of stationary of the total potential.

The variational principle exists whether the operator specifying the mixed form is

symmetric or self-adjoint [see 61, Cap.11]. In this case not only the displacements

but also the stresses and/or strains are used as primary variables. The mixed finite

elements existing in literature come from the combination of the various extended

variational principles with the use of many different finite element interpolation

schemes [see 19, Section.4.4].

While displacement-based finite elements are in general most popular, there are

some situations, such as the analysis of incompressible media and the analysis of

plate and shell structures, where the use of mixed finite is clearly more effective [see

19, pg. 276].

Analogously, concerning the modeling of limit states (collapse or shakedown),

many different approaches are reported with kinematic, static or mixed finite ele-

ments being used. This diversity of finite elements comes from the combination of

any of the principles presented in Sect. 2.3 with different proposals of interpolation.

The definition of the “best”strategy of interpolation applicable to all shakedown

problems is an extremely difficult task and, in our opinion, it seems to remain an

open issue, at least in the near future.

In this work structural finite elements are developed based on a mixed interpola-

tion scheme with stress and velocity fields independently approximated by standard

Lagrange polynomial functions. Hence, for elastic problems the proposed elements

correspond to the family of finite elements obtained from the Hellinger-Reissner

functional. Although there is no claim that this interpolation strategy provides the

best results in terms of accuracy and/or performance, our choice was supported by

the following reasons:

21

(a) The excellent performance demonstrated by this type of mixed finite element

for the solution of the shakedown problems in plane strain conditions [10, 60];

(b) The discretization of the kinematic principle (2.41) with kinematic elements

result in unbounded dissipation D(d), unless the interpolated velocity fields v

only produce isochoric strain rates [10, 21]. This restriction in general leads to

the construction of either laborious or very limited interpolation fields.

(c) It is reported that the exact fulfillment of equilibrium constraints in the statical

formulation (2.39) is hard to achieve with static elements, except for special loads

[21]. Moreover, statical elements are built using non-standard finite element

procedures making their implementation cumbersome.

It was shown in Chapter 2 that the ideally elastic stress fields σE of the structure

are important input data in shakedown formulation. Approximations for σE can

be obtained by means of linear analyses. In this case we employ the same finite

elements and mesh used in the calculation of shakedown to compute the elastic

stresses for each extreme loading assuming that the structure is made of an ideally

elastic material. In theory this is not the unique way for obtaining the elastic

stresses, however, this approach is very convenient specially because there is no

effort in post-processing elastic analysis results for the preparation of the input

data for shakedown analysis. More importantly, the use of the same finite elements

and mesh for obtaining the elastic stresses and for the solution of the shakedown

problem makes sure that the approximated elastic stresses which are input data

for the shakedown analysis are fully compatible with the interpolation functions

employed to solve the shakedown problem. As a consequence, the equilibrium of

elastic stresses with the amplified external loads during the shakedown analysis is

fulfilled.

Accuracy of approximated elastic stress fields has impact on the shakedown so-

lution. In particular to the cases where the failure mechanism is alternate plasticity,

the elastic stress estimate plays a critical role due to the local nature of this phe-

nomenon. Hence, in this Chapter the focus is on the properties that mixed finite

elements should present to provide good approximations for the elastic problem.

We start in the continuum setting showing the variational principle and the

regularity required for the trial functions. Then, we present the general strategy of

discretization used in the three finite elements proposed in this work for the elastic

problem. We remark that details are provided in the respective chapters of each

finite element.

Finally, some remarks on stability and convergence in elasticity for mixed finite

elements are briefly outlined, mainly following [61].

22

3.1 Two field mixed formulation in elasticity

Consider the Hellinger-Reissner functional [20, 61]

ΠHR(σ, u) = −1

2

∫Bσ · E−1σdB +

∫Bσ · DudB −

∫Bu · bdB −

∫Γτ

u · τ dΓτ . (3.1)

then, invoking the stationarity of its first variation (δΠHR(σ, u) = 0) one can arrive

at the following variational or weak formulation for linear elastic problems:

Given the prescribed body forces b acting on B and the boundary tractions τ

applied on the part Γτ of the boundary Γ of B, find the displacement field u ∈ Vand stress field σ ∈ W ′ that satisfy∫

Bδσ · DudB −

∫Bδσ · E−1σdB = 0 ∀ δσ ∈ W ′ (3.2)

and ∫Bσ · DδudB −

∫Bδu · bdB −

∫Γτ

δu · τ dΓτ = 0 ∀ δu ∈ V (3.3)

The space V consists of all vector fields sufficiently regular and satisfying homo-

geneous constraints in the part Γu, of the boundary Γ of B. Note that since the first

derivative of u and its variation δu appears in (3.2) and (3.3), respectively, elements

of V should have integrable derivatives, i.e., u, δu ∈ C0. On the other hand, the

field spaceW ′ comprises less regular functions. In fact, the elements ofW ′ are sym-

metric tensor fields which can be piecewise continuous functions, i.e., σ, δσ ∈ C−1.

The interpolation functions must comply with the above requirements.

3.2 Discrete formulation in elasticity

3.2.1 Interpolation functions

For interpolation schemes compatible with the Hellinger-Reissner functional, displa-

cement and stress fields can be given by the following general form:

u(x) ≈ u(x) = Nvui, (3.4)

σ(x) ≈ σ(x) = Nσσi (3.5)

where: Nv(x) and Nσ(x) denote the interpolation functions for displacement and

stress respectively, and ui and σi are the vectors which collect the displacement and

23

stress nodal coefficients. Analogously, for the fields of variations we have:

δu(x) ≈ δu(x) = Nvδui, (3.6)

δσ(x) ≈ δσ(x) = Nσδσi. (3.7)

The strain rate field is calculated for each point x in the body by using

ε(x) ≈ ε(x) = DNvui (3.8)

where the symbol D, in this case, denotes the discrete tangent deformation operator.

3.2.2 Discrete equilibrium, compatibility and constitutive

equations

In the sequel we derive the discrete equilibrium, compatibility and constitutive equa-

tions starting from the Hellinger-Reissner functional.

Firstly, the fields in (3.2)-(3.3) are replaced by their respective counterparts

defined in (3.4)-(3.7) leading to(∫BNTσDNvdB

)ui · δσi −

(∫BNTσ E−1NσdB

)σi · δσi = 0 ∀δσi, (3.9)

and(∫BNTσDNvdB

)δvi · σi −

(∫BNTv bdB +

∫Γτ

NTv τ dΓτ

)· δvi = 0 ∀δvi. (3.10)

Then, since δσi and δvi are arbitrary vectors, equations (3.9) and (3.10) can be

recast in the following matrix form

Bui − E−1σi = 0 (3.11)

and

BTσi − F = 0. (3.12)

by defining the matrices

F :=

∫BNTv bdB +

∫Γτ

NTv τ dΓτ , (3.13)

E−1 :=

∫BNTσ E−1Nσ dB, (3.14)

B :=

∫BNTσDNvdB, (3.15)

24

where the symbol F denotes the load vector, the symbol E corresponds to the discrete

elastic relation; and the symbol B represents the discrete strain operator which is

the same operator defined in (2.38) for the shakedown problem. These matrices are

calculated and assembled element by element. Their final forms are slightly different

for the three finite elements presented in this work. Details of each implementation

are provided in the three chapters describing each finite element.

Note that (3.12) is the discrete equilibrium equation and (3.11) is the discrete

constitutive equation. This latter is rearranged as

σi = EBui = Eεi. (3.16)

The discrete compatibility equation between displacement and strain fields is

εi = Bui (3.17)

provided the definition of the vector εi ∈ Rq containing the nodal strain rate com-

ponents.

Finally, by combining (3.12) and (3.16) one can obtain the classical matrix re-

presentation

Kui = F (3.18)

where K is the elastic stiffness matrix given by

K := BTEB. (3.19)

3.2.3 Remarks on stability and convergence in elasticity

In addition to the regularity required for the interpolation functions discussed in

Sect. 3.1, mixed u − σ finite elements must comply with the following important

requirement

nσ ≥ nu (3.20)

where nσ denotes the number of unknowns in terms of stress components and nu

refers to the number of unknowns in displacements (number of degrees of freedom).

The above constraint is called the count condition and it is a necessary although

not sufficient condition to ensure stability in linear elastic problems [61]. Locking

or non-convergent results are expected when the count condition is violated.

To obtain good results in linear elastic analysis by using mixed u − σ finite

elements care must be taken with some other issues, such as: the fulfillment of

patch test requirements, the enforcement of excessive (meaningless) continuity on

interpolated stress fields, the observance of the frame invariance requirement (results

25

independent of the orientation of the initially chosen coordinate system), etc. To a

deeper discussion on this matter we refer to ZIENKIEWICZ et al. [see 61, Cap.11].

26

Chapter 4

Solution Algorithm for Shakedown

Analysis

The set of optimality conditions (2.48)-(2.53) presented in Chapter 2 is solved using

a tailor-made algorithm, which was originally developed for limit analysis problems

considering non-linear yield functions [60, 66], and was further improved for sha-

kedown analysis. This algorithm is briefly outlined in what follows because details

are not essential to the objectives of this work. Indeed, the finite elements proposed

here are fully applicable to other solution techniques. Therefore, it is our intention

to keep the exposition of the algorithm as compact as possible. For details we refer

to [10, 22].

4.1 The linearization process

Equations from [10, 22] are simplified for materials obeying the von Mises yield

criterion. Internal variables or dependence on temperature are disregarded. Then

the yield function takes the form

f(σ) :=3

2‖ S ‖2 (4.1)

where S denotes the deviatoric stress tensor defined as

S := σdev = σ − σm1 = Pσ (4.2)

with

σm :=1

3trσ =

1

31 · σ (4.3)

P :=I− 1

31 ⊗ 1 (4.4)

27

The symbols σm and P respectively denote the mean stress and the fourth order

tensor projecting stress tensors over the subspace of deviators. The symbols 1 and

I are the second and fourth order identity tensors.

As a consequence the gradient and Hessian of the yield function used in the

algorithm are

∇σf(σ) = 3S, ∇σσf(σ) = 3P . (4.5)

Optimality conditions (2.48)-(2.53) are written in matrix form by means of the

following vectors and matrices

λ ∈ Rm, λ :=λk; k=1:m

(4.6)

j ∈ Rm, jk := σk · ∇σf(σF + µσk + σr) (4.7)

F ∈ Rm×m, F := diag(f(σF + µσk + σr)

)(4.8)

G ∈ Rq×m, colkG := ∇σf(σF + µσk + σr) (4.9)

where λ is the vector of plastic multiplier components, F has zero entries except

in its diagonal which entries are given byf(σF + µσk + σr); k=1:m

, and colkG

identifies the column k of matrix G.

With the above definitions the discrete optimality conditions are recast in the

compact vector equation

g(x) :=

BTσr

Bv − G(µ,σr)λ

1− jG(µ,σr) · λ−FG(µ,σr)λ

=

0

0

0

0

(4.10)

f(σF + µσk + σr) ≥ 0, λk ≥ 0 k=1:m. (4.11)

where the vector of unknowns is x :=[v σr µ λk

]T.

The vector equation (4.10) is solved for a givenσk; k=1:m

and σF under

the restrictions given by (4.11). Matrices j(µ,σr), F(µ,σr) and G(µ,σr) are now

identified as non-linear matrix functions of the unknown variables.

A Newton-type iterative process is employed to solve (4.10) with each iteration

computing the transformation

x→ xnew = x + δx (4.12)

where x denotes the present value, xnew the next iterate and δx the increment.

28

Then, the explicit increment is

v →vnew, (4.13)

σr →σrnew = σr + δσr, (4.14)

µ→µnew = µ+ δµ, (4.15)

λ→λ0. (4.16)

The Newton iteration is −∇g(x)(xnew − x) = g(x) with the gradient of g(x)

given by

−∇g(x) =

0 −BT 0 0

−B H b G0 bT c jT

0 ΛG Λj F

(4.17)

with

H :=∑

λk∇σσf(σF + µσk + σr) ∈ Rq×q, (4.18)

b :=∑

λk∇σσf(σF + µσk + σr)σk ∈ Rq, (4.19)

c :=∑[

∇σσf(σF + µσk + σr)σk · σk]∈ R, (4.20)

Λ :=diag(λk) ∈ Rm×m. (4.21)

The Newton-like prediction is made by solving the Newton iteration for the

unknownsvnew, δσr, δµ, λ

0

applying the step increment as defined above.

4.2 The solution algorithm

The iteration process of the solution algorithm consists of two basic steps: (i) a

Newton iteration for the set of equalities (4.10); (ii) the enforcement of plastic

admissibility, given by the inequalities (4.11), by means of relaxation and uniform

scaling of the variables.

4.2.1 Newton iterations

The Newton iteration process is performed taking into account the favorable struc-

ture of matrices obtained with the finite element method. This procedure given

in detail in [10] leads to significant reduction in computational cost. Matrices are

computed element by element, as usually made in the finite element method.

The original linear system is transformed by a sequence of clever substitutions

resulting in the following steps [10]:

29

1. Compute, element by element:

(a) the matrix E ∈ Rq×q of fictitious elastic coefficients

colkG =∇σf(σF + µσk + σr) k=1:m (4.22)

jk =σk · ∇σf(σF + µσk + σr) k=1:m (4.23)

H =∑

λk∇σσf(σF + µσk + σr) (4.24)

b =∑

λk∇σσf(σF + µσk + σr)σk (4.25)

c =∑[

∇σσf(σF + µσk + σr)σk · σk]

(4.26)

−Λ−1F =diag

(f(σF + µσk + σr)

λk

)(4.27)

Q =H−1G (4.28)

W =GTH−1G −Λ−1F (4.29)

E =H−1 −QW−1QT (4.30)

(b) the fictitious stress

σ = QW−1j (4.31)

(c) the contribution of each element to the pseudo-energy

c = c = W−1j · j − Eb · b− 2σ · b (4.32)

(d) the contribution of each element to the stiffness and modifying loadings

K =BT EB (4.33)

F µ =BT(Eb+ σ

)(4.34)

F r =−BTσr (4.35)

2. Solve the global linear system

Kvµ =F µ (4.36)

Kvr =F r (4.37)

30

3. Update the vector of unknowns

δµ =1− F µ · vr

c+ F µ · vµ(4.38)

vnew =δµvµ + vr (4.39)

λ0

=W−1QT (Bvnew − δµb) + δµW−1j (4.40)

δσr =H−1(Bvnew − δµb− Gλ0

)(4.41)

4.2.2 Relaxation and uniform scaling of the variables

An important feature of the solution procedure is that it consist of an interior point

algorithm. To this end plastic admissibility constraints (4.11) are imposed at the

end of each iteration.

Plastic admissibility is enforced by means of a factor p and a step relaxation

factor s such that

σnew = pY s, Y s := (σr + sδσr) (4.42)

µnew = pµs, µs := (µ+ sδµ) . (4.43)

For the determination of the scaling factor p the maximum amplification that

still satisfies plastic admissibility at the end of the step is computed for each is plastic

admissibility constraint assuming an initial guess for the step relaxation factor s.

The stress state at the end of an iteration is strictly admissible if

f(σF + µnewσk + σr

new)< 0. (4.44)

Thus

f(σF + p(µsσ

k + Y s))≤ γff(σF + µσk + σr) < 0 (4.45)

where γf is a parameter defined as

γf = min

γ0f ,δµ

µ

(4.46)

and γ0f is a control parameter arbitrarily chosen in the interval (0, 1). Then, mini-

mum of the scaling factors for all elements is used.

The procedure for the determination of the scaling factor p must be repeated

if the new estimate µnew for the shakedown safety factor is lower than the present

estimate µ [10]. In this case a smaller relaxation factor s taken from a prescribed

sequence such as 1, s, s2, ... is used.

After the determination of the scaling factor p and the relaxation factor s the

31

residual stress vector σr and the safety factor µ are updated. Note that the velocity

vector is updated by vnew = δµvµ+vr. The last parameter that must be updated is

the vector of plastic multipliers λ. This operation is performed tanking into account

strict positiveness of λ resulting in an invertible matrix Λ. The applied rule is(λnew

)k← max

(λ0)k, γλ ‖ λ0 ‖∞

(4.47)

where γλ is a prescribed tolerance and ‖ · ‖∞ is the norm of the maximum absolute

value of components.

32

Chapter 5

Continuum-based Beam Element

In this chapter the continuum-based (CB) approach [12, 13, 67] is used for the

development of a three-node beam element with rectangular cross-section for the

solution of 2-D shakedown problems.

This is the first finite element developed during this research with the aim of

testing the combination of a mixed interpolation scheme with the continuum-based

approach for collapse and shakedown calculation by direct methods. The methodo-

logy tested for this simpler 2-D beam was lately employed in the development of the

axisymmetric shell element presented in Chapter 6 and finally in the formulation of

the 3-D shell element described in Chapter 7.

The beam element employs a interpolation strategy with displacements and ve-

locities quadratic along the element axis and linear through its thickness. In each

layer, stresses are interpolated linearly with respect to mid-surface parameters and

also linearly across the layer thickness. Continuity of stresses through layers may or

may not be enforced, giving two alternative versions of the finite element.

In the following sections we provide details for the implementation of this CB

beam element and the numerical examples presented in [68].

5.1 Geometry, displacements and strains

The global (fixed) orthonormal reference frame is denoted R = ex, ey. We also

use ez = ex ∧ ey to represent plane rotations.

33

Figure 5.1: Continuum-based beam element

5.1.1 Master and slave nodes

In the one-dimensional element i there are three master nodes and three correspon-

ding directors, respectively denoted as:

xa = xaex + yaey, a=1:3 (5.1)

pa = cos θaex + sin θaey, a=1:3 (5.2)

The underlying continuum element i has six slave nodes:

xα = xαex + yαey, α=s1:s6. (5.3)

For convenience, we use the following two alternative label sequences in order to

identify the slave nodes of element i (see Figure 5.1):

s1:s6 ≡ 1−, 1+, 2−, 2+, 3−, 3+. (5.4)

Using this notation we can relate master and slave nodes by the following formulas

xa− = xa − ha

2pa, a=1:3, (5.5)

xa+ = xa + ha

2pa, a=1:3, (5.6)

where ha is the beam thickness at node a.

The fiber through the master node a and parallel to the director pa (a pseudo

normal) moves rigidly. Thus, the velocities va, va− and va+ of master and slave

34

nodes are related by the angular velocity ωaez of the director pa as follows

va− = va − ha

2ωaez ∧ pa, a=1:3, (5.7)

va+ = va + ha

2ωaez ∧ pa, a=1:3. (5.8)

Considering (5.6) and (5.5) the above leads to

va− = va + ωaez ∧(xa− − xa

), (5.9)

va+ = va + ωaez ∧(xa+ − xa

). (5.10)

We gather (5.9) and (5.10) in symbolic matrix notation as[va−

va+

]= T i,a

[va

ωa

]where T a,i :=

[1 2 ez ∧ (xa− − xa)

1 2 ez ∧ (xa+ − xa)

](5.11)

with 1 2 denoting the 2-D identity tensor.

Relations (5.11) are now collected, for a=1:3, in one intrinsic equation

vs,i = T ivi (5.12)

by defining the block-diagonal matrix

T i := diag(T 1,i,T 2,i,T 3,i

)(5.13)

and

vs,i :=[vs1 vs2 vs3 vs4 vs5 vs6

]i,T(5.14)

≡[v1− v1+ v2− v2+ v3− v3+

]i,T. (5.15)

We recall that the superscript T is used to denote the transpose of a matrix.

For the sake of completeness, we write now the intrinsic equation (5.12) in global

coordinates. In the present notation this reads as follows

[vs,i]R

=[T i]R [

vi]R, (5.16)

35

with

[vs,i]R

=[vs1x vs1y . . . vs6x vs6y

]R,T(5.17)

=[v1−x v1−

y v1+x v1+

y . . . v3−x v3−

y v3+x v3+

y

]R,T, (5.18)

[vi]R

=[v1x v1

y ω1 . . . v3x v3

y ω3]R,T

, (5.19)

and

[T a,i

]R=

1 0 ya − ya−

0 1 xa− − xa

1 0 ya − ya+

0 1 xa+ − xa

R

. (5.20)

5.1.2 Kinematics of the underlying continuum element

The parent element (see Figure 5.1) is

Ω := (ξ, η) ∈ [−1, 1]× [−1, 1], (5.21)

and the geometry mapping in the i-th finite element domain Bi is

Ω 3 ξ = (ξ, η) 7→ x(ξ)|Bi =∑

α=s1:s6

gα(ξ)xα (5.22)

with the symbol gα(ξ) denoting the Lagrange interpolation functions

g1 = 14ξ(ξ − 1)(1− η), (5.23)

g2 = 14ξ(ξ − 1)(1 + η), (5.24)

g3 = 12(1− ξ)(1 + ξ)(1− η), (5.25)

g4 = 12(1− ξ)(1 + ξ)(1 + η), (5.26)

g5 = 14ξ(ξ + 1)(1− η), (5.27)

g6 = 14ξ(ξ + 1)(1 + η). (5.28)

Note that (5.22) is intrinsic.

The curvilinear coordinate system R(ξ) = ex, ey, used to enforce the hypothe-

sis of zero transverse normal stress is obtained employing the same method used by

BELYTSCHKO et al. [12]. Accordingly, we define the base vectors ex tangent and

36

ey normal to the lamina (see Figure 5.1) as

ex :=x,ξ‖x,ξ‖

=x,ξe

x + y,ξey√

(x,ξ)2 + (y,ξ)2, (5.29)

ey := ez ∧ ex =−y,ξex + x,ξe

y√(x,ξ)2 + (y,ξ)2

, (5.30)

with

x,ξ(ξ) =∑

α=s1:s6

gα,ξxα, y,ξ(ξ) =

∑α=s1:s6

gα,ξyα. (5.31)

The subscript “, ξ”denotes the corresponding derivative.

Change of basis from R(ξ) to R is accomplished by using the general formulas

for vectors and second order tensors

[a]R = R [a]R or [A]R = R [A]R RT (5.32)

with

R(ξ) :=

[ex · ex ex · ey

ey · ex ey · ey

]. (5.33)

The Jacobian of the geometry mapping is given by

J(ξ) = ∇Tx(ξ) =∑

α=s1:s6

xα ⊗∇gα(ξ) (5.34)

or [x,ξ x,η

y,ξ y,η

]R=

∑α=s1:s6

[xαgα,ξ xαgα,η

yαgα,ξ yαgα,η

]R. (5.35)

Further, the chain rule gives

gα,ξ = gα,xx,ξ + gα,yy,ξ (5.36)

gα,η = gα,xx,η + gα,yy,η (5.37)

or

∇ξgα = JT∇xgα, (5.38)

and thus

∇xgα = J−T∇ξgα. (5.39)

Finally, the derivatives of the interpolation functions with respect to spatial

coordinates are obtained by substituting the inverse of the Jacobian in (5.39), i.e.[gα,x

gα,y

]R=

1

x,ξy,η − x,ηy,ξ

[y,η −y,ξ−x,η x,ξ

]R [gα,ξ

gα,η

](5.40)

37

for α = s1 : s6.

5.1.3 The interpolation of displacements and velocities

The assumed displacement and velocity fields are then

u(ξ)|Bi =∑

α=s1:s6

gα(ξ)uα, (5.41)

v(ξ)|Bi =∑

α=s1:s6

gα(ξ)vα. (5.42)

Since gα(ξ)uα = [gα(ξ)1 2]uα, we can recast the above interpolations in the

following compact intrinsic form

u(ξ)|Bi = Nv(ξ)us,i, v(ξ)|Bi = Nv(ξ)vs,i, (5.43)

where the vector vs,i, collecting velocities of all slave nodes, was defined in (5.14)

and

us,i =[us1 us2 us3 us4 us5 us6

]i,T(5.44)

≡[u1− u1+ u2− u2+ u3− u3+

]i,T, (5.45)

Nv(ξ) :=[g11 2 g21 2 g31 2 g41 2 g51 2 g61 2

]. (5.46)

Note that (5.43), (5.44), (5.14) and (5.46) are intrinsic and that the component

representations of all vectors, u(ξ), v(ξ), us,i and vs,i, depend on the adopted

coordinate system. For instance, we wrote the global representation of vs,i in (5.17).

However, the linear map Nv is represented by the same matrix in any reference

system

[Nv(ξ)]R =

[g1 0 g2 0 g3 0 g4 0 g5 0 g6 0

0 g1 0 g2 0 g3 0 g4 0 g5 0 g6

]. (5.47)

We can now express the velocity field, written in curvilinear components, i.e.

[v(ξ)]R, in terms of the parameters of interpolation, which in this case are the nodal

velocities of the element written in global components, i.e. [vi]R

. To this end we

first write (5.43) in curvilinear components, i.e.

[v(ξ)]R(ξ) = [Nv(ξ)][vs,i]R(ξ)

. (5.48)

38

Next, we define the matrix

R(ξ) := diag(R,R,R,R,R,R) (5.49)

with R(ξ) given in (5.33).

Therefore, according to (5.32), the following holds

[vs,i]R(ξ)

= RT (ξ)[vs,i]R. (5.50)

Further, by introducing (5.12) in the above,

[vs,i]R(ξ)

= RT (ξ)[T i]R [

vi]R. (5.51)

Finally, by combining (5.48) and (5.51) we get

[v(ξ)]R(ξ) = [Nv(ξ)]RT (ξ)[T i]R [

vi]R. (5.52)

5.1.4 Enforcing bending theory hypotheses

We assume that in the direction perpendicular to the local laminar axis both strain

and stress components are zero. Accordingly, we use the following notation for

planar tensors

d = dxex ⊗ ex + dye

y ⊗ ey + d(xy)

[1√2

(ex ⊗ ey + ey ⊗ ex

)](5.53)

and

σ = σxex ⊗ ex + σye

y ⊗ ey + σ(xy)

[1√2


)]. (5.54)

Thus, we have that

d(xy) =√

2dxy, σ(xy) =√

2σxy, (5.55)

and the standard Euclidean structure applies, i.e.

σ · d = σxdx + σydy + σ(xy)d(xy), (5.56)

‖d‖2 = d2x + d2

y + d2(xy), (5.57)

‖σ‖2 = σ2x + σ2

y + σ2(xy). (5.58)

39

The relevant compatibility equations are then

dx = vx,x, (5.59)

d(xy) = 1√2

(vx,y + vy,x) . (5.60)

5.1.5 Computing strain in a generic point

We consider next the computation of the infinitesimal strain rate ∇sym v at a generic

(Gauss) point of the underlying continuum. Remarks by BELYTSCHKO et al. [12,

p.104] are relevant to this matter.

Let x(ξk) denote a generic (Gauss) point. According to (5.51)

[vs,i]R(ξk)

= RT (ξk)[T i]R [

vi]R

(5.61)

is the column vector of velocity components, of all slave nodes, in the curvilinear

system of the generic point considered. Moreover, we write, using (5.48)

[v(ξ)]R(ξk) = [Nv(ξ)][vs,i]R(ξk)

. (5.62)

This equation is simply the componentwise counterpart of the intrinsic interpolation

formula (5.42) when the fixed global reference system is chosen coincident with the

curvilinear frame R(ξk).

Consequently, we can adopt the following approach in order to compute the

strain rate tensor in the generic point x(ξk), written in curvilinear components: (i)

we use (5.61) to compute curvilinear nodal velocity components, in R(ξk), from the

general interpolation parameters (the cartesian nodal velocity components) and (ii)

we use (5.62) and the standard procedure for strain computations in isoparametric

finite elements. That is,

1. Find the curvilinear coordinate system R(ξk) = ex, ey given by (5.29) and

(5.30) and the rotation matrix R(ξk) using (5.33).

2. Compute, by using (5.61), all slave node velocities in the curvilinear system

of the generic point

[vs,i]R(ξk)

= RT (ξk)[T i]R [

vi]R.

3. For α=s1:s6 apply the change of coordinates for slave node positions

xα = RT (ξk) xα. (5.63)

4. For α=s1:s6 find (gα,x, gα,y).

40

To this end, compute the Jacobian in the curvilinear frame using (5.35) with

(xa, ya) (from (5.63)) instead of (xa, ya). Then, use (5.39) written for curvili-

near coordinates, i.e. analogous to (5.40).

5. Compute the infinitesimal strain tensor

[∇sym

x v(ξk)]R(ξk)

=[∇sym

x Nv(ξk)]R(ξk) [

vs,i]R(ξk)

. (5.64)

Here (5.64) is written, discarding zero strain components, as

[vx,x

1√2

(vx,y + vy,x)

]R(ξk)

=[G1 . . .G6

]R(ξk)

vs1x

vs1y...

vs6x

vs6y

R(ξk)

(5.65)

with auxiliary matrices Gα defined by

Gα(ξ) :=

gα,x 0

1√2gα,y

1√2gα,x

R , α=s1:s6. (5.66)

In summary, the above procedure explains how to compute the compact formula

[∇symx v(ξ)]R(ξ) = [∇sym

x Nv(ξ)]R(ξ)RT (ξ)[T i]R [

vi]R, (5.67)

giving the strain rate, in local laminar components, at a generic point x(ξ) and in

terms of the interpolation parameters [vi]R.

5.2 Stresses

Layers are determined by a partition of the parent domain Ω =⋃`=1:L Ω` given as

follows: First, we select a sequence of coordinates −1 = η0, η1, . . . , ηL = 1. Then,

each parent layer is set as Ω` = (−1, 1)× (η`−1, η`). Spatial layers are the image of

parent layers in the geometry mapping (5.22).

Let us define now a local coordinate η ∈ (−1, 1), to be used in the restricted

domain of one single layer `, given by

η :=η − a`b`

, (5.68)

where a` := 12

(η` + η`−1) and b` := 12

(η` − η`−1).

41

We choose to interpolate the bending normal stress and shear components in

each layer ` using bilinear functions. Thus, each layer has four stress nodes located

at corners and

σx(ξ)|Bi =∑j=1:4

tj(ξ)σj/`x , (5.69)

σ(xy)(ξ)|Bi =∑j=1:4

tj(ξ)σj/`(xy), (5.70)

where the interpolation parameters σj|`x and σ

j|`(xy) are interpreted as the stress com-

ponents of the j-th stress node of layer ` in its own laminar directions. The inter-

polation functions are

tj(ξ) = 14(1 + ξjξ)(1 + ηj η) (5.71)

or

tj(ξ) = 14(1 + ξjξ)

[1 +

ηj − a`b2`

(η − a`)], (5.72)

with (ξj, ηj) denoting the coordinates of node j. These are the usual bilinear func-

tions, just restricted to one layer, where an appropriate transversal coordinate was

defined in (5.68). This change of variable is also used to compute integrals by

numerical quadrature.

The above interpolation of stress is cast in matrix form as

[σ(ξ)]R = [Nσ(ξ)][σ`]

(5.73)

with (1 2 is the identity matrix 2× 2)

[Nσ(ξ)] = [t11 2 . . . t41 2] (5.74)

by defining in each layer `

[σ(ξ)]R =[σx σ(xy)

]T, (5.75)[

σ`]

=[σ1|` σ2|` σ3|` σ4|`]T , (5.76)

and in each stress node t of layer `

[σt|`]

=[σt|`x σ

t|`(xy)

]T. (5.77)

42

5.3 Yield function

Consider the intrinsic yield function associated to the von Mises yield criterion for

a continuum body (excluding dependence on temperature or internal variables)

f(σ) =3

2‖ S ‖ −σ2

Y ≥ 0 (5.78)

where the symbol S denotes the deviator stress tensor

S := σ − 1

3trσ1 2. (5.79)

The equation 5.78 above is rewritten in laminar components as

f(σ) =1

2

[(σx − σy)2 + (σy − σz)2 + (σz − σx)2

+ 3(σ2(xy) + σ2

(yz) + σ2(zx))

]− σ2

Y ≥ 0 (5.80)

The yield function for the proposed beam element is easily obtained by elimina-

ting from (5.80) all the null stress components. Thus,

f(σ) = σ2x +

3

2csσ

2(xy) − σ2

Y , (5.81)

The parameter cs is added to the model to consider or not the contribution of the

shear stresses to yielding. Accordingly, this coefficient can be set equal to 0 or 1. If

cs = 0 the classical beam theory is simulated, otherwise, cs = 1 and shear stresses

are considered.

The gradient and the Hessian of the yield function are then

∇σf =

[2σx

3csσ(xy)

]and ∇σσf =

[2 0

0 3cs

](5.82)

5.4 Discrete strain operator

We now proceed to define the discrete strain operator. To this end we consider the

work-conjugacy relation for an arbitrary element i:

δW int =

∫Bi

[σ]R · [∇symx δv]R rdx. (5.83)

In the discrete setting we may define for each layer ` a work-conjugate strain

rate vector [d`]

:=[d1|` d2|` d3|` d4|`

]T, (5.84)

43

where [dt|`]

=[dt|`x d

t|`(xy)

]T(5.85)

and t denotes the stress nodes of layer `.

With the above definitions, (5.83) can be re-written in discrete form as

δW int =∑`

[σ`]·[δd`], (5.86)

where[δd`]

is given by [δd`]

= B`,i[δvi]R

(5.87)

and B`,i denotes the layer discrete strain operator obtained as follows.

In terms of parent element coordinates (5.83) is expressed as

δW int =

∫Ω

[σ(ξ)]R · [∇symx δv(ξ)]R detJ dξ dη, (5.88)

and can be re-written as sum of element layer contributions

δW int =∑`

∫Ω`

[σ(ξ)]R · [∇symx δv(ξ)]R detJ dξ dη. (5.89)

We now use (5.67) and (5.73) in (5.88). Then, after some manipulations we get

δW int =∑`

[σ`]·∫Ω`

[Nσ(ξ)]T [∇symx Nv(ξ)]R(ξ)RT (ξ) detJ dξ dη

[T i]R [

δvi]R

(5.90)

From the above, the layer discrete strain operator is identified as

B`,i :=

∫Ω`

[Nσ(ξ)]T [∇symx Nv(ξ)]R(ξ)RT (ξ) detJ dξ dη

[T i]R. (5.91)

We recall that the nodal strain rate components defined in (5.84) and (5.85) are

approximations of the corresponding strain multiplied by the layer volume.

Further, let us define

A`,i(ξ) := [Nσ(ξ)]T [∇symx Nv(ξ)]R(ξ) . (5.92)

44

Then, considering (5.66) and (5.73), we have

A`,i(ξ) =

t11 2

t21 2

t31 2

t41 2

[G1 G2 . . . G5 G6

]R. (5.93)

Thus

A`,i(ξ) =[As1,`,i . . .As6,`,i

]R, (5.94)

with

Aα,`,i =

t11 2

t21 2

t31 2

t41 2

[Gα

]R, (5.95)

or, explicitly,

Aα,`,i =

t1gα,x 0

1√2t1gα,y

1√2t1gα,x

......

t4gα,x 0

1√2t4gα,y

1√2t4gα,x

R

. (5.96)

Next, we define

B`,islave :=

∫Ω`

[As1,`,i . . .As6,`,i

]RRT detJ dξ dη

(5.97)

and substitute (5.97) in (5.91) to compute the strain operator for the layer ` of the

i-th element as

B`,i := B`,islave

[T i]R. (5.98)

Finally, the element strain operator is obtained by assembling the contributions

of each layer:

Bi =∑`

B`,i. (5.99)

The assembly in (5.99) is easily performed by means of well-known procedures

employed in finite element analysis and depends on whether or not continuity of

stress fields between layer is enforced. When stress fields are discontinuous between

45

layers, the assembly is simply

Bi =

B1,islave...

BL,islave

[T i]R. (5.100)

5.5 Discrete elastic relation

This section is just concerned with obtaining the discrete form, denoted E, of the

elastic operator for to the CB beam element proposed for shakedown analysis. It

is worth remembering here that this matrix is only used in the data preparation

for the analysis, when the discrete approximation of the ideally elastic stress fields

σE, for each extreme loading, are computed. This elastic pre-analysis gives rise to

the practical definition of the prescribed domain of loading variations introduced in

(2.20) and (2.32). Clearly, obtaining the discrete operator E can be performed using

the Hellinger-Reissner principle (see (3.14)). However, we adopt a shorter equivalent

presentation described in the following.

Consider the elastic response of the beam, where both ideally-elastic stresses and

strains are related by

σE(x) = E εE(x) ∀x ∈ B. (5.101)

The elastic stress-strain operator E above is defined for the CB beam considering

the plane stress σy = 0 constraint. Then,

[E]R

= E

[1 0

0 11+ν

], (5.102)

where E and ν are Young’s modulus and Poisson’s ratio.

The elastic relationship (5.101) is substituted, for the i-th element, by the con-

dition ∫Biδσ · εE dx =

∫Biδσ · E−1 σE dx (5.103)

that holds for all stress variations δσ compatible with the adopted stress interpola-

tion. Also εE and σE are solutions of the approximate formulation, so that∫Biδσ · E−1 σE dx =

∑`

∫Ω`

[δσ(ξ)]R ·[E−1

]R [σE(ξ)

]R(ξ) detJ dξ dη, (5.104)

and substituting (5.73) in the previous equation we conclude that∫Biδσ · E−1 σE dx =

∑`

[δσ`]·(E`,i)−1 [

σ`], (5.105)

46

where (E`,i)−1

:=

∫Ω`

[Nσ(ξ)]T[E−1

]R[Nσ(ξ)] detJ dξ dη (5.106)

is the inverse of the elastic relation for the layer ` of the i-th element.

By performing analogous operations on the right hand side of (5.103), we arrive

at ∑`

[δσ`]·B`,i

[ui]R − (E`,i)−1 [

σ`]

= 0, ∀δσ`. (5.107)

Then, for each element layer, the following relation holds

B`,i[ui]R

=(E`,i)−1 [

σ`], (5.108)

and, adopting the symbol [σi] to denote the stress components of the i-th element,

it is also true that

Bi[ui]R

=(Ei)−1 [

σi]

(5.109)

with

Bi =∑`

B`,i and(Ei)−1

=∑`

(E`,i)−1

, (5.110)

where the symbols Bi and Ei respectively represent for the i-th element the discrete

strain operator and the discrete elastic matrix.

Finally, the discrete stress-strain relation for the i-th element, can be written as

[σi]

= EiBi[ui]R. (5.111)

The discrete elastic relation Ei, in turn, can be generally calculated in two steps:

firstly, its inverse (Ei)−1 is obtained by assembling the contributions (5.106) of each

element layer ` using standard assembly procedures. Secondly, the assembled matrix

is inverted. Note, however, that the computational burden to obtain the matrix Ei is

greatly reduced when there is no imposed continuity on stress fields between layers.

In this case, the matrices (E`,i)−1 of each layer are calculated and inverted and then

assembled directly in the matrix Ei.

5.6 Numerical Examples

Three numerical examples are provided for the proposed CB beam element with

different purposes. The objective of the first example is the assessment of the element

accuracy in linear elastic analysis of thick and thin beams. The aim of the second

example is the verification of the beam element in the computation of collapse loads.

In the last example the beam element is verified in shakedown analysis.

Analytical solutions are used in the three examples as reference to results obtai-

47

Figure 5.2: Semi-circular arch subjected to a load 2P

ned with the CB beam element.

5.6.1 Linear elastic solution of a curved beam

The first problem, depicted in Fig. 5.2, consists of a semi-circular arch subjected to

a vertical load 2P . The beam has a rectangular cross-section with width b, height

h, area Ab = bh and moment of inertia I = bh3/12. Displacements and rotation at

points A and B (see Fig. 5.2) calculated by using Castigliano’s energy theorem are

uany (A) = −(

3π − 8

4

)PR3

EI− πPR

4EAb− πPR

4GAbScc, (5.112)

uanx (B) =PR3

2EI− PR

2EAb+

PR

2GAbScc, (5.113)

ωan(B) =PR2

EI

(π2− 1). (5.114)

The applied load is P = 1.0. The material and geometric properties are set to

the following values: Young Modulus E = 1.0 × 109, Poisson’s ratio ν = 0.3, shear

modulus G = E/[2(1 + ν)], shear correction coefficient for a rectangular section

Scc = 0.85, radius R = 10 and width b = 1. Four slenderness ratios R/h are tested:

10, 100, 700 and 1000. Due to symmetry, only a half of the structure needs to be

modeled with proper boundary conditions. A mesh refinement study is conducted

by applying four finite element meshes with increasing degree of discretization: 1, 2,

5 and 10 elements. The elements have only one layer for stress interpolation. The

third term on the RHS of (5.112) and (5.113) are the part of the displacement due

to shear. These terms are usually negligible for thin beams.

Figure 5.3 summarize the results for the four slenderness ratios and for the four

meshes. The graph on the left presents the non-dimensional displacement uy/uany at

the center of the arch (point A in Fig. 5.2). The graphs on the center and on the right

48

exhibit respectively the non-dimensional displacement ux/uanx and rotation ω/ωan at

the support (point B in Fig. 5.2). Note that accuracy increases with discretization.

The displacement error remains almost stable when slenderness ratio increases and

no locking is observed. Meshes with more than two elements present error less than

0.8%.

If shear stresses are neglected then the analytical expressions for tangential

uant (φ) and normal uann (φ) displacements of the arch centerline are:

uant (φ) =1

2R[2a · sinφ+ (2− π)a · cosφ− (a+ b)φ cosφ− 2aφ− (2− π)a] (5.115)

uann (φ) =1

2R[(a+ b) cosφ− (2− π)a · sinφ+ (a+ b)φ sinφ− 2a] (5.116)

where, a =R2P

EIand b =

P

EAb(5.117)

The above equations are used in Fig. 5.4 to verify the deformed shape of the arch for

the finite element models with R/h = 700. Three distinct FE meshes with increasing

degree of discretization are considered. It can be seen a very good agreement between

analytical and numerical results.

5.6.2 Collapse of a beam in tension plus bending

As an example of limit analysis consider a beam with rectangular cross section (h

is the height and b is the width) under uniform external axial load N and external

bending moment M , as sketched in Fig. 5.5.

The axial stress σx is a scalar function of the transversal coordinate z and the

equilibrium equations are:

N = b

∫ h/2

−h/2σxdz, M = −b

∫ h/2

−h/2σxzdz (5.118)

The material of the beam is perfectly plastic with yield stress σY . Plastic ad-

missibility is determined by:

P : −σY ≤ σx(z) ≤ σY (5.119)

The collapse limits for pure traction and bending are, respectively, NY = bhσY

and MY = bh2σY /4.

Considering that the system of loads is defined by the following non-dimensional

load parameters n = N/NY and m = M/MY and recalling that the loading (nc,mc)

producing instantaneous plastic collapse satisfies n2c + mc = 1, hence, the collapse

49

Figure 5.3: Nodal displacements and rotation for the semi-circular arch consideringvarious slenderness ratios

50

Figure 5.4: Deformation of the circular arch: analytical solution and FE results

Figure 5.5: Beam under axial traction and bending moment

51

Figure 5.6: Collapse limits for a beam section under axial traction and bendingmoment: analytical solution and FE results. Left: discontinuous stress field betweenlayers. Right: continuous stress field between layers.

factor α for the load (n,m) is determined analytically by solving

(αn)2 + αm = 1 (5.120)

Figure 5.6 depicts a comparison between the analytical solution given by (5.120)

and the results obtained with the CB beam element. The structure is modeled

with just one element and different number of layers: 2, 4, 8 and 16. The number

of layers are increased to assess the effect of enriching interpolated stress fields on

element results. The graph on the left presents the results obtained assuming that

the stress field is discontinuous between adjacent layers whereas the graph on the

right depicts the results obtained when continuity of stresses is enforced. It can be

seen that, the higher the number of layers the closer numerical and analytical results

are, as expected. For elements with 8 or more layers numerical and analytical results

are quite close.

Considering a given number of layers, elements with no continuity of stress fields

between layers are in general more accurate. For the elements with 8 layers the

maximum difference to the analytical collapse loads is less than 1% for the element

with discontinuous stress fields and less than 3% for the element with continuous

stress fields (see Tab. 5.1). This better performance of elements with no continuity

of stress is also expected since axial stress distribution through the thickness is not

continuous in the collapse of a elastic-perfectly plastic beam subjected to bending

plus traction.

52

Table 5.1: Maximum differences between analytical and FE resultsNumber of Elem. discont. stress field Elem. cont. stress field

Layers No. stress nodes Max. error [%] No.stress nodes Max. error [%]2 8 11.1 6 33.34 16 2.8 10 8.48 32 0.7 18 2.116 64 0.3 34 0.5

On the other hand, elements without continuity of stress field between adjacent

layers are computationally much more expensive than the respective version with

enforced continuity between layers. Indeed, when the number of layers increases

the number of stress nodes of an element with discontinuous stress fields between

layers tends to be twice the number of stress nodes of an respective element with

continuous stress fields between layers. The stress nodes are the points where plastic

admissibility is verified by the optimization algorithm for shakedown analysis. Then,

if the number of stress nodes increases the computational cost increases as well.

It is worth emphasizing that for this example the accuracy of results of elements

with similar number of stress nodes are comparable as shown in Tab. 5.1.

5.6.3 Shakedown analysis in tension plus bending

The structure considered in this example is exactly the same presented in the previ-

ous example except for the loads which are slightly modified to allow a shakedown

analysis. In this case, the axial load N remains fixed whereas the external bending

moment M can vary, i.e., −M ≤ M ≤ M . The solution for this problem is well

known [10, 69] and is provided here in the form of the Bree diagram shown in Figure

5.7. The blue curve labeled with C in Figure 5.7 indicates the collapse limit for the

loads already given in the previous example. The red curve represents the elastic

shakedown limits, and the green line is the pure elasticity boundary. Labels AP

and IC indicate respectively alternating plasticity and incremental collapse, which

are the failure mechanisms for loads beyond elastic shakedown limits. The objective

of this example is to compare the analytical limit for elastic shakedown with finite

element results obtained with the mixed CB beam element.

The numerical models for this example are identical to those employed for limit

analysis. However, in this case only three number of layers are considered: 2, 4 and

8. The finite element results are presented in Figure 5.8. The graph on the left in

Figure 5.8 depicts results for the element with discontinuous stresses between layers.

The graph on the right, in turn, shows the results for the element with continuous

stress field within the element. Note that all simulations present the same outcome.

Indeed, even the element with 2 layers and continuity of stress, which was poor

53

Figure 5.7: Bree diagram for a beam under constant axial traction and a variableuniform bending moment, ZOUAIN [10].

for the determination of accurate collapse limit, is capable of capturing the correct

limits for shakedown in this particular case.

54

Figure 5.8: Bree diagram for shakedown analysis of a beam under constant axialtraction and a variable uniform bending moment:analytical solution and FE re-sults. Left: discontinuous stress field between layers. Right: continuous stress fieldbetween layers.

55

Chapter 6

Continuum-based Axisymmetric

Shell Element

In this chapter a three-node shell element is devised for the solution of axisymmetric

shakedown problems.

This element is analogous to the beam element presented in the previous chapter.

Indeed, due to the mathematical similarity between the axisymmentric problem and

the plane stress (or plane strain) problem, the equations for this axisymmetric shell

element can be readily derived from the equations of the CB beam element by

considering that [see 61, Cap.5]:

1. In a body with symmetry of revolution in geometry and loads, the two com-

ponents of displacement in any plane section passing through its axis of sym-

metry define completely the state of strain and, consequently, the stresses on

the body;

2. The global reference frame of an axisymmetric body is given by the radial and

axial coordinates of the point. By using this reference frame instead of the

global reference frame employed for the CB beam element, one can realize that

the interpolation functions of the axisymmetric element and beam element are

identical;

3. All integrations of the axisymmetric element must consider the volume defined

by rotating the area of the element around its axis of revolution;

4. In axisymmetric bodies any displacement in radial direction induces strains

and stresses in circumferential direction. Therefore, a third component of

strain and stress in circumferential direction must be considered in the axisym-

metric shell element. This is the most remarkable difference between the

axisymmetric finite element and the CB beam element.

56

Figure 6.1: Continuum-based axisymmetric shell element: displacement nodes (mas-ter and slave) and stress nodes

Although there are in literature various works focusing on the limit analysis of

axisymmetric structures, references on shakedown analysis are rare and shell ele-

ments are seldom used. Some numerical examples of shakedown analysis of axisym-

metric structures can be found, for instance, in [41, 70–76], most of them employing

only solid elements. In this work the outcomes of the proposed axisymmetric shell

element [77] are compared to the results from solid finite element models and avai-

lable analytical solutions after the detailed presentation of the element formulation.

In what follows we organize the sections with the same divisions used in the

previous chapter. Some intermediate steps of calculation presented for the CB beam

element are omitted in this chapter to make the text more concise. However, it is

our intention to make the implementation of each finite element easier by giving

in its respective chapter the equations in their final form. For this reason and for

the sake of clarity of the exposition, some formulae previously introduced may be

repeated in this chapter.


The global (fixed) orthonormal reference frame (see Fig 6.1) is denoted R =

er, ez ≡ ex1 , ex2 with ez corresponding to the axis of symmetry. The circum-

ferential direction, perpendicular to the axisymmetric half-plane, is represented by

eϕ = er ∧ ez.

We remark that equations describing the geometry and the displacements of this

CB axisymmetric element are exactly the same expressions presented in Sect. 5.1.1–

5.1.3 to define the geometry and the displacements of the CB beam element with

57

coordinates x and y replaced by r and z respectively.


In the one-dimensional element i there are three master nodes and three correspon-

ding directors, respectively denoted as:

xa = raer + zaez, a=1:3, (6.1)

pa = cos θaer + sin θaez, a=1:3. (6.2)

The respective underlying continuum element has six slave nodes:

xα = rαer + zαez, α=s1:s6. (6.3)

Using the notation introduced in (5.4) we can also relate master and slave nodes

by the formulas (5.5) and (5.6), i.e.,

xa− = xa − ha

2pa, a=1:3, (6.4)

xa+ = xa + ha

2pa, a=1:3, (6.5)

where ha is the shell thickness at node a.

A fiber through the master node a and parallel to the director pa (a pseudo

normal) moves rigidly. Thus, the velocities va, va− and va+ of master and slave

nodes are related by the angular velocity ωaeϕ of the director pa through the same

equations (5.7)-(5.10) used in the CB beam formulation

va− = va − ha

2ωaeϕ ∧ pa ≡ va + ωaeϕ ∧

(xa− − xa

), a=1:3, (6.6)

va+ = va + ha

2ωaeϕ ∧ pa ≡ va + ωaeϕ ∧

(xa+ − xa

), a=1:3. (6.7)

or in the symbolic notation (5.11)[va−

va+

]= T i,a

[va

ωa

]where T a,i :=

[1 2 eϕ ∧ (xa− − xa)

1 2 eϕ ∧ (xa+ − xa)

]. (6.8)

Relations between velocities of master and slave nodes can be also collected in

the intrinsic equation (5.12) which is given in global coordinates as

[vs,i]R

=[T i]R [

vi]R, (6.9)

58

with

[vs,i]R

=[vs1r vs1z . . . vs6r vs6z

]R,T(6.10)

=[v1−r v1−

z . . . v3+r v3+

z

]R,T, (6.11)

[vi]R

=[v1r v1

z ω1 . . . v3r v3

z ω3]R,T

, (6.12)

and

[T a,i

]R=

1 0 za − za−

0 1 ra− − ra

1 0 za − za+

0 1 ra+ − ra

R

. (6.13)


The geometry of the CB axissymetric shell element is mapped with the same equa-

tions employed for the CB beam element. The parent element domain is defined by

Ω := (ξ, η) ∈ [−1, 1]× [−1, 1] and the geometry mapping in the i-th finite element

domain Bi is given by (5.22),i.e.,

Ω 3 ξ = (ξ, η) 7→ x(ξ)|Bi =∑

α=s1:s6

gα(ξ)xα, (6.14)

with the symbol gα(ξ) denoting the Lagrange interpolation functions g1 = (1/4)ξ(ξ−1)(1− η), g2 = (1/4)ξ(ξ − 1)(1 + η), g3 = (1/2)(1− ξ)(1 + ξ)(1− η), g4 = (1/2)(1−ξ)(1 + ξ)(1 + η), g5 = (1/4)ξ(ξ + 1)(1− η) and g6 = (1/4)ξ(ξ + 1)(1 + η).

The curvilinear coordinate system R(ξ) = ex, ey, used to enforce the hy-

pothesis of zero transverse normal stress, is obtained by defining the base vectors

ex tangent and ey normal to the lamina (see Figure 6.1) as

ex :=x,ξ‖x,ξ‖

=r,ξe

r + z,ξez(

r2,ξ + z2

,ξ

)1/2, (6.15)

ey := eϕ ∧ ex =−z,ξer + r,ξe

z(r2,ξ + z2

,ξ

)1/2, (6.16)

where subscript “, ξ”denotes the corresponding derivative.


for vectors and second order tensors (5.32) with

R(ξ) :=

[er · ex er · ey

ez · ex ez · ey

]. (6.17)

59

The Jacobian of the geometry mapping (J(ξ) in (5.34)) is given by[r,ξ r,η

z,ξ z,η

]R=

∑α=s1:s12

[rαgα,ξ rαgα,η

zαgα,ξ zαgα,η

]R. (6.18)

and the derivatives of the interpolation functions with respect to spatial coordinates

are obtained by substituting the inverse of the Jacobian in the equation (5.39), i.e.[gα,r

gα,z

]R=

1

r,ξz,η − r,ηz,ξ

[z,η −z,ξ−r,η r,ξ

]R [gα,ξ

gα,η

](6.19)

for α=s1:s6.


The assumed displacement and velocity fields are interpolated using exactly the same

equations employed for the CB beam. The expressions presented in Sect. 5.1.3 are

not repeated here because they can be used for the CB axisymmetric shell element

without any modification.


We assume that in the direction perpendicular to the local laminar axis both strain

and stress components are zero and keep all other non-zero strain components pos-

sible in an axisymmetric deformation. Accordingly, we use the notation for planar

tensors introduced in Sect. 5.1.4 to write

d = dxex ⊗ ex + dϕe

ϕ ⊗ eϕ + d(xy)

[1√2


)](6.20)

and

σ = σxex ⊗ ex + σϕe

ϕ ⊗ eϕ + σ(xy)

[1√2


)]. (6.21)

The relevant compatibility equations for the CB axisymmetric shell element are

dx = vx,x, (6.22)

dϕ =vrr, (6.23)

d(xy) = 1√2

(vx,y + vy,x) . (6.24)

Note that the component dϕ does not appear in the CB beam element formulation.

The velocity in (global) radial direction vr in (6.23) is computed from the velo-

60

city components in curvilinear coordinates applying (5.32) with Ri,j given in (6.17).

Accordingly

vr = R1,1vx + R1,2vy. (6.25)

6.1.5 Computing strain in a generic point

The computation of the infinitesimal strain rate at a generic (Gauss) point of the un-

derlying continuum and in terms of the interpolation parameters [vi]R is performed

with the formula (5.67)



vi]R.

Note that in this equation the vector of strain rates ∇symx v is given in local la-

minar components whereas the vector of nodal velocities [vi]R is written in global

components.

The above equation is calculated by following the same procedure explained for

the CB beam element in subsection 5.1.5, i.e.,

1. Find the curvilinear coordinate system R(ξk) = ex, ey given by (6.15) and

(6.16) and the rotation matrix R(ξk) using (6.17);

2. Compute, by using

[vs,i]R(ξk)

= RT (ξk)[T i]R [

vi]R

all slave node velocities in the curvilinear system of the generic point;


xα = RT (ξk) xα; (6.26)

4. For α=s1:s6 find (gα,x, gα,y);


(xa, ya) (from (6.26)) instead of (ra, za). Then, use (5.39) written for curvili-

near coordinates, i.e. analogous to (6.19).

5. Finally, compute the infinitesimal strain tensor

[∇sym

x v(ξk)]R(ξk)

=[∇sym

x Nv(ξk)]R(ξk) [

vs,i]R(ξk)

. (6.27)

61

Here (6.27) is written, discarding zero strain components, as

vx,xvrr

1√2

(vx,y + vy,x)

R(ξk)

=[G1 . . .G12

]R(ξk)

vs1x

vs1y...

vs6x

vs6y

R(ξk)

(6.28)


Gα(ξ) :=

gα,x 0

R1,1

rgα

R1,2

rgα

1√2gα,y

1√2gα,x

R

, α=s1:s6. (6.29)

6.2 Stresses

Layers are determined by a partition of the parent domain Ω =⋃`=1:L Ω` given in the

same way as in the CB beam element, i.e.: first, we select a sequence of coordinates

−1 = η0, η1, . . . , ηL = 1. Then, each parent layer is set as Ω` = (−1, 1)× (η`−1, η`).

Spatial layers are the image of these L parent layers through the geometry mapping

(6.14).

A local coordinate η ∈ (−1, 1), to be used in the restricted domain of one single

layer `, is defined as

η :=η − a`b`

, (6.30)

where a` := 12

(η` + η`−1) and b` := 12

(η` − η`−1).

Bending normal stress and shear components in each layer ` are interpolated

with bilinear functions. Each layer has four stress nodes located at corners and

σx(ξ)|Bi =∑j=1:4

tj(ξ)σj/`x , (6.31)

σϕ(ξ)|Bi =∑j=1:4

tj(ξ)σj/`ϕ , (6.32)

σ(xy)(ξ)|Bi =∑j=1:4

tj(ξ)σj/`(xy), (6.33)

where the interpolation parameters σj|`x , σ

j|`ϕ and σ

j|`(xy) are interpreted as the stress

components of the j-th stress node of layer ` in its own laminar directions. Note that

the component σj|`ϕ , not present in the CB beam element, is a remarkable difference

in the CB axisymmetric shell formulation.

62

The interpolation functions are the same employed in the CB beam, then

tj(ξ) = 14(1 + ξjξ)(1 + ηj η) ≡ 1

4(1 + ξjξ)

[1 +

ηj − a`b2`

(η − a`)], (6.34)

with (ξj, ηj) denoting the coordinates of node j. These are the usual bilinear func-

tions, just restricted to one layer, where an appropriate transversal coordinate was

defined in (6.30). This change of variable is also used to compute integrals by

numerical quadrature.

The above interpolation of stress is cast in matrix form as

[σ(ξ)]R = [Nσ(ξ)][σ`]

(6.35)


[Nσ(ξ)] = [t11 3 . . . t41 3] (6.36)


[σ(ξ)]R =[σx σϕ σ(xy)

]T, (6.37)[

σ`]

=[σ1|` σ2|` σ3|` σ4|`]T , (6.38)


[σt|`]

=[σt|`x σt|`ϕ σ

t|`(xy)

]T. (6.39)

6.3 Yield function

The yield function f(σ) for the axisymmetric CB shell element is obtained from

the von Mises yield criterion (5.78) by considering the non-zero stress components.

Thus,

f(σ) = σ2x + σ2

ϕ − σxσϕ +3

2csσ

2(xy) − σ2

Y . (6.40)

The parameter cs is added to the model to consider or not the contribution of the

shear stresses to yielding. This coefficient which is set equal 1 by default can also

be set equal to 0 to disregard the shear stress component.


∇σf =

2σx − 1

2σϕ − 1

3csσ(xy)

and ∇σσf =

2 −1 0

−1 2 0

0 0 3cs

(6.41)

63


For obtaining the discrete strain operator for the axisymmetric CB shell element

we firstly define in the discrete setting for each layer ` a work-conjugate strain rate

vector [d`]

:=[d1|` d2|` d3|` d4|`

]T, (6.42)

where [dt|`]

=[dt|`x dt|`ϕ d

t|`(xy)

]T(6.43)


Then, we follow the steps presented for the CB beam element in Sect. 5.4 to

identify the layer discrete strain operator as

B`,i := B`,islave

[T i]R

(6.44)

where

B`,islave :=

2π

∫Ω`

[As1,`,i . . .As6,`,i

]RRT r(ξ) detJ dξ dη

(6.45)

with the auxiliary matrix A defined as


Then, by considering (6.29) and (7.57), the matrix A is computed as

A`,i(ξ) =

t11 3

t21 3

t31 3

t41 3

[G1 G2 . . . G5 G6

]R. (6.47)

Alternatively,

A`,i(ξ) =[As1,`,i . . .As6,`,i

]R, (6.48)

with

Aα,`,i =

t11 3

t21 3

t31 3

t41 3

[Gα

]R, (6.49)

64

or, explicitly,

Aα,`,i =

t1gα,x 0

R1,1

rt1gα

R1,2

rt1gα

1√2t1gα,y

1√2t1gα,x

......

t4gα,x 0

R1,1

rt4gα

R1,2

rt4gα

1√2t4gα,y

1√2t4gα,x

R

. (6.50)


of each layer:

Bi =∑`

B`,i. (6.51)



stress fields between layers is enforced. When stress fields are discontinuous between


Bi =

B1,islave...

BL,islave

[T i]R. (6.52)


In this section we compute the discrete form of the elastic operator for the axisymme-

tric CB shell element denoted by E. This matrix is only used in the data preparation




(2.20) and (2.32).

Firstly, we define the elastic stress-strain operator E considering both plane stress

σz = 0 and axisymmetric constraints. Then,

[E]R

=E

1− ν2

1 ν 0

ν 1 0

0 0 1− ν2

, (6.53)


Hence, by following the demonstration steps presented for the CB beam element

65

Figure 6.2: Cylindrical shell under ring load

in section 5.4, we conclude that the inverse of the discrete elastic operator for an

element i is computed by assembling the contributions of each layer ` as

(Ei)−1

=∑`

(E`,i)−1

, (6.54)

where(E`,i)−1

denotes the inverse of the elastic relation for the layer ` of the i-th

element defined as

(E`,i)−1

:= 2π

∫Ω`

[Nσ(ξ)]T[E−1

]R[Nσ(ξ)] r(ξ) detJ dξ dη (6.55)

The discrete elastic relation Ei, can be generally calculated in two steps: firstly,

its inverse (Ei)−1 is obtained by assembling the contributions (6.55) of each element

layer ` using standard assembly procedures. Secondly, the assembled matrix is

inverted. Note, however, that the computational burden to obtain the matrix Ei is

greatly reduced when there is no imposed continuity on stress fields between layers.

In this case, the matrices (E`,i)−1 of each layer are calculated and inverted and then

assembled directly in the matrix Ei.


6.6.1 Limit analysis of cylindrical shells under a ring load

In this example we consider cylindrical shells, with different lengths, under the action

of a ring load applied in their central cross-section. The model is geometrically

defined by the semi-length of the cylinder L, the mean radius R and the shell

thickness h, as shown in Figure 6.2. The material of the shell is assumed to be

elastic-ideally plastic with von Mises yield criterion and yield stress σY .

66

Analytical solutions for this example were obtained by various authors, e.g. DE-

MIR [78], SAWCZUK and HODGE [79], ZOUAIN [59] and CHAKRABARTY [80],

by assuming different simplified yield functions in generalized variables. In the fol-

lowing, we discuss the qualitative differences of these approaches with respect to the

present CB shell modeling.

It was generally identified in the technical literature that these ring-loaded cy-

linders can be classified in long, medium and short tubes according to the following

types of collapse mechanism, exhibited under plastic collapse conditions:

Long cylinders are those whose collapse mechanism presents a central circum-

ferential plastic hinge surrounded by two symmetric zones undergoing distributed

plastic strain rates and two rigid zones at both ends. Two additional symmetric

plastic hinges, separating the central plastic zone and the rigid endings, are found

in analytical solutions based on uniform or sandwich shell assumptions combined

with Mises or Tresca plastic conditions. However, these symmetric plastic hinges

are not present in our CB shell solutions; a fact confirmed in our additional solution

using solid rotationally symmetric finite elements.

Medium cylinders exhibit one central plastic hinge and distributed plastic defor-

mation all along their length. No rigid regions are observed.

For short cylinders there is no plastic hinge and plastic collapse happens with

dissipation distributed along the whole tube.

In the analytical solution given by [78] the cylinder was idealized as a uniform

shell made of elastic-ideally plastic material with Tresca yield criterion. The mecha-

nisms of collapse were identified in [78], depending on the non-dimensional length

of the cylinder

ξL =L√Rh

. (6.56)

Accordingly, for uniform Tresca shells, ring-loaded cylinders can be classified as

follows:

ξL ≥ 2.336↔ long cylinder, (6.57)

0.725 < ξL < 2.336↔ medium cylinder, (6.58)

ξL ≤ 0.725↔ short cylinder. (6.59)

The above categories are used here for the definition of the geometries of the solid

models employed to validate the proposed CB shell element.

Earlier analytical solutions for the particular case of long cylinders are found in

[79]. These solutions were obtained by the combination of two different shell models

– the uniform shell and the sandwich shell – with two different yield criteria – the

von Mises and the Tresca criterion. Lately, ZOUAIN [59] provided the analytical

67

solution for cylinders of different lengths considering a sandwich shell model and

the Tresca yield criterion. More recently, this problem of cylinder under ring load

was revisited in [80] where analytical solutions were given for sandwich and uniform

shells considering different yield criteria. A recent application of this example in the

verification of a numerical method for limit analysis is found in [81].

We remark that all analytical solutions aforementioned were obtained under the

assumption of Kirchhoff-Love shell theory where fibers initially normal to the mid-

plane of the shell are regarded to remain plane and normal after deformation. This

hypothesis, generally used for thin shells, differs from the kinematic assumptions of

the the proposed CB axisymmetric shell element (see Sect. 6.1.4).

In the present example, finite element approximations are obtained for cylin-

ders with radius to thickness ratio R/h = 30. The ring load is applied as a non-

dimensional reference load f0, same as that used in [78], i.e.,

f0 =

√M0

RN0

=1

2

√h

R, (6.60)

where

N0 = hσY and M0 =h2

4σY (6.61)

are, respectively, the limit loads corresponding to axial force and pure bending.

Due to symmetry, only a half of the cylinder is modeled with n-layer CB axisym-

metric shell elements without continuity of stresses between layers. Mesh density as

well as the number of layers are determined by a mesh sensitivity study comprising

two steps:

(i) a coarser mesh made of 4-layer elements is refined by adding more elements

until convergence;

(ii) the number of layers of the finer mesh defined in the previous step is gradually

increased until a new convergence.

The accuracy of CB shell results is verified by solid models built with 6-node

mixed triangular axisymmetric elements with their limit loads calculated using the

same algorithm employed with the CB shell models [10, 22]. In this verification,

cylinders with four different lengths are considered: a short (ξL = 0.511), a transition

short-medium (ξL = 0.741), a medium (ξL = 1.716) and a long (ξL = 3.536) cylinder.

The interpolation of the triangular solid element is made in a similar fashion

than the CB shell element, i.e. displacements are assumed quadratic and continuous

whereas stresses are linear and discontinuous. This mixed interpolation scheme has

been successfully used in shakedown analyses as well as other applications [22, 82].

68

Figure 6.3: Finite element meshes

The mesh density of solid models is also determined by a mesh refinement proce-

dure to guarantee convergence. As a result, regions of the model presenting higher

plastic strain rates are meshed with a higher density of elements. In Fig. 6.3

is illustratively given the coarser and the finer meshes used for the long cylinder

(ξL = 3.536) model. Figure 6.3a shows the meshes for the CB shell and Fig. 6.3b

depicts the meshes employed in solid models.

To illustrate the convergence pattern obtained in both mesh refinement proce-

dures, it is shown in Fig. 6.4 two graphs corresponding to the long cylinder case.

Figure 6.4a depicts the limit loads computed in step (i) of the mesh refinement pro-

cedure for the CB shell, whereas Fig. 6.4b shows the results obtained in step (ii). In

both graphs the results of the respective solid models with different levels of mesh

refinement are also presented for comparison.

Next, we show in Fig. 6.5 the different collapse mechanisms computed with the

CB shell element for the long, medium and short cylinders.

Figure 6.6 depicts an example of plastic strain rate plot, which corresponds to

instantaneous collapse of the short cylinder (ξL = 0.511); there we can see the plastic

strain rates distributed along the length of the cylinder, as expected.

The load multipliers (µ = F/f0), calculated with both shell and solid models,

are provided in Table 6.1. Note that both results are very close, with maximum

difference of 1%. However, the difference in number of degrees of freedom between

shell and solid models makes clear the advantage of using the CB shell element in

69

Figure 6.4: Example of mesh refinement results: long cylinder (ξL = 3.536)

Figure 6.5: Collapse mechanisms of cylinders under ring load

70

Figure 6.6: Plastic strain rates for a short cylinder (ξL = 0.511) under collapse dueto a ring load

this case.

Table 6.1: Load multipliers (µ) obtained with CB shell models and solid axisymme-tric models.

CB shell models Solid models

ξL stress stress

d.o.f. nodes µ d.o.f. nodes µ

0.511 49 192 1.000 1,581 1,104 0.999

0.741 91 360 1.337 8,613 6,240 1.341

1.716 127 504 1.684 7,477 5,376 1.701

3.536 127 672 1.899 15,093 11,004 1.918

The collapse load computed analytically by SAWCZUK and HODGE [79] for

a long cylinder, considering a uniform shell and the von Mises yield criterion, was

µ = 1.949. This result is 2.63% higher than the collapse load calculated in the

present work with the CB shell. However, the analytical solution was computed in

[79] neglecting the effect of shear stresses. Thus, that result can be considered as

less precise than the collapse loads obtained here, with both shell and solid models.

Moreover, by setting to zero the shear coefficient cs in the yield function of the CB

elements, (5.81), we obtain a result comparable to the solution from [79]. In fact,

with cs = 0 the collapse load predicted by the CB shell model is µ = 1.938 which is

only 0.56% less than the analytical result from [79].

Finally, we plot in Fig. 6.7 the collapse load as a function of the cylinder length.

This graph is constructed employing the CB shell element with 8 layers. The figure

also shows, for comparison purpose, the results of the four solid models and also

the analytical solutions given in [78] and [59], the latter calculated considering the

Tresca’s yield criterion and respectively uniform or sandwich shells.

71

Figure 6.7: Collapse load versus non-dimensional cylinder’s length

6.6.2 Limit and shakedown analysis of pressure vessels with

ellipsoidal and torispherical heads

In this example we employ the CB shell element to calculate collapse and shakedown

limits for pressure vessels with ellipsoidal and torispherical heads under internal

pressure. The geometry and material properties are taken from [72, 76]. In both

references the material was assumed as elastic-perfectly plastic with von Mises yield

criterion and axisymmetric solid elements were used with different techniques to

calculate the limit loads.

The pressure vessels have the same thickness in the cylindrical part and in the

head. Their geometry, schematically given in Fig. 6.8, is defined by non-dimensional

parameters. For pressure vessels with ellipsoidal heads these parameters are a/b and

a/h. For pressure vessels with torispherical heads the parameter r/a is also used.

The radius R of the spherical part of the torispherical head can be calculated

with the recursive equation

R = b+√

(r + c−R) (r − c−R) (6.62)

where the parameter c := a−r denotes the distance from the the axis of the cylinder

to center of the torus. Note that R > r + c.

The transition sphere-torus occurs at the critical radius given by

rcrit = c

[1 +

(R

r− 1

)−1]. (6.63)

72

Figure 6.8: Geometrical parameters of pressure vessels

The geometries considered in [76] were:

• for ellipsoidal head

a/b = 2, 2.5, 3, 4

a/h = 10, 25, 100, 300

• for torispherical head

a/b= 1.0, 1.5625, 2.0, 2.5, 3.33333

r/a= 0.12, 0.20, 0.30

a/h= 10, 100/7, 25, 50

and the the finite element models were carried out in [76] with the ABAQUS code

neglecting large displacement effects. The material properties employed were the

Young modulus E = 70 GPa, the Poisson’s ratio ν = 0.3 and the yield stress

σY = 20 MPa.

The estimate of the collapse pressure, sometimes cumbersome when calculated

via incremental analyses, was made in [76] by monitoring the pressure-deflection

73

curve. The collapse pressure was defined in [76] as the applied pressure from which

a very small increase in pressure produced a large increase in displacement.

The computation of shakedown loads, even more difficult by incremental analy-

sis, was made in [76] by using an iterative procedure. The methodology used by

[76] consisted of varying the load between zero and an arbitrary pressure p and mo-

nitoring in each cycle the displacement of the top of the head, the plastic strain,

and the residual stress field. If the displacement was repeated at the next cycle and

the plastic strain did not change during a cycle, and the residual stress field at zero

pressure did not violate the yield condition, then the pressure vessel was regarded as

shaking down and the pressure p was increased, otherwise, the pressure was reduced.

The procedure was repeated until the difference between pressures to shake down

and not to shake down showed less than 1.0% according to that authors.

The influence of the size of knuckle radius r in failure mechanism of a torispherical

head was investigated in [72]. The parameter R/a was used instead of parameter

a/b. Two cases were examined: r/a = 0.1 and r/a = 0.4, both with R/a = 1.5 and

a/h = 25.

Instead of using incremental analyses, it was employed in [72] a direct approach

to calculate the limit loads. The method used in [72] consisted of finding the collapse

and shakedown loads making use of Melan’s static theorem. The formulation was

posed as a conic quadratic optimization problem. The pressure vessel heads were

discretized in [72] with 6-node triangular axisymmetric solid elements using uniform

meshes.

The shakedown analyses of [72] demonstrated that alternate plasticity was the

critical failure mode for the pressure vessel with lower knuckle radius while the other

pressure vessel was more prone to fail due to incremental or instantaneous collapse.

This conclusion corroborates the fact that the knuckles with small radii are elastic

stress concentrators that can lead to alternate plasticity failure of the pressure vessel.

No material properties were listed in [72]. Nevertheless, for this example in

particular, the relevant material property is the Poisson’s ratio, assumed here as

ν = 0.3. The exact Young modulus is not important in this case because the

hypothetically elastic (unbounded) stress field only depends on the Poisson’s ratio,

the boundary conditions and the applied mechanical loads. The exact value of yield

stress is also irrelevant because the limit pressure is linearly dependent on the yield

stress and the applied internal pressure p is normalized with respect to the reference

value

p0 = σYh

a. (6.64)

In this work, the finite element models are carried out with the CB shell elements

with 8 layers and stress fields discontinuous between adjacent layers. Mesh density

74

as well as number of layers are defined using a mesh refinement procedure equivalent

to that described in the previous example (see Sect. 6.6.1).

Limit loads calculated with the CB shell models are listed in Table 6.2 for pres-

sure vessels with elliptical heads and in Table 6.3 for pressure vessels with torisphe-

rical heads.

Table 6.2: Limit loads calculated with CB shell element for pressure vessels withellipsoidal heads. pE = elastic limit; pLA = collapse load; pSD = shakedown load;Mch = failure mechanism; AP = alternate plasticity; IC = incremental collapse.

a/b a/h pE pLA pSD Mch

2 10 0.679 1.181 1.181 IC2 25 0.673 1.134 1.131 IC2 100 0.661 1.053 1.053 IC2 300 0.657 0.942 0.941 IC

2.5 10 0.482 1.034 0.964 AP2.5 25 0.454 0.910 0.908 AP2.5 100 0.415 0.708 0.707 IC2.5 300 0.394 0.584 0.583 IC

3 10 0.369 0.878 0.739 AP3 25 0.341 0.731 0.682 AP3 100 0.296 0.527 0.526 IC3 300 0.270 0.417 0.416 IC

4 10 0.185 0.694 0.370 AP4 25 0.229 0.527 0.457 AP4 100 0.187 0.355 0.355 IC4 300 0.162 0.266 0.266 IC

Besides the models of Table 6.2, we also perform some additional limit analyses

for of pressure vessels with ellipsoidal heads. In these analyses we assume the ge-

ometrical parameter a/b = 2 fixed and consider very high a/h ratios (= 1500, 10

000 and 100 000). The objective of such models is just to verify the CB shell ele-

ment response when its thickness approaches zero. The results presented in Fig. 6.9

compare very well with the results from [76], also calculated using axisymmetric

shell elements from ABAQUS. No locking was observed and both results approach

asymptotically to the collapse limit pLA = 0.7559 based on membrane theory [83]

as the thickness decreases.

Limit analysis of the set of pressure vessels with ellipsoidal heads shows different

mechanisms of instantaneous collapse (see Fig. 6.10). In general, the failure of these

structures occurs in the head. The exception is the pressure vessel with a/b = 2

and a/h = 10 with the failure taking place in the head and in the cylindrical part

simultaneously. For moderately thick pressure vessels, 10 ≤ a/h ≤ 25, the plastic

strain rates spread throughout the head with higher values taking place in the head

vertex. The majority of the pressure vessels with very thin wall thicknesses a/h ≥

75

Table 6.3: Limit loads calculated with CB shell element for pressure vessels withtorispherical heads. pE = elastic limit; pLA = collapse load; pSD = shakedown load;Mch = failure mechanism; AP = alternate plasticity; IC = incremental collapse.

a/b r/a a/h pE pLA pSD Mch

1 1 25 1.120 1.163 1.163 IC1.563 0.2 25 0.631 1.164 1.164 IC1.563 0.3 14.29 0.763 1.175 1.175 IC

2 0.12 25 0.401 1.098 0.803 AP2 0.2 25 0.447 1.092 0.894 AP2 0.3 14.29 0.539 1.148 1.077 AP2 0.3 25 0.505 1.073 1.01 AP

2.5 0.12 10 0.361 1.108 0.723 AP2.5 0.12 25 0.305 0.840 0.610 AP2.5 0.12 50 0.252 0.629 0.503 AP2.5 0.2 14.29 0.372 0.948 0.743 AP2.5 0.2 25 0.331 0.797 0.662 AP2.5 0.2 50 0.289 0.628 0.578 AP2.5 0.3 10 0.378 0.803 0.756 AP2.5 0.3 14.29 0.344 0.729 0.687 AP2.5 0.3 25 0.306 0.646 0.612 AP2.5 0.3 50 0.268 0.544 0.537 AP

3.333 0.12 10 0.262 0.803 0.523 AP3.333 0.12 14.29 0.248 0.722 0.495 AP3.333 0.12 25 0.217 0.583 0.433 AP3.333 0.12 50 0.178 0.433 0.356 AP3.333 0.2 14.29 0.236 0.565 0.472 AP3.333 0.2 25 0.200 0.474 0.400 AP3.333 0.2 50 0.167 0.372 0.335 AP

R/a r/a a/h pE pLA pSD Mch

1.5 0.1 25 0.319 0.921 0.638 AP1.5 0.4 25 0.718 1.164 1.164 IC

100 presented the failure with plastic strain rates distributed in the intersection of

the head and the cylindrical part.

The mechanisms of instantaneous collapse found for the set of pressure vessels

with torispherical head are similar to those described for the pressure vessels with

ellipsoidal head. Depending on the geometry, the weakest part can be the head or

the cylinder. In this case, the radius of the knuckle has an important effect on the

head strength. To illustrate, we provide in Fig. 6.11 the collapse mechanisms of the

models investigated in [72]. It can be seen that for the pressure vessel with r/a = 0.1

the collapse occurs in the head whereas for the pressure vessel with r/a = 0.4 the

cylinder is the weakest part.

The difference between the CB shell results and the results from literature [72, 76]

is presented in Table 6.4 for the pressure vessels with ellipsoidal head and in Table

76

Figure 6.9: Collapse load versus thickness ratio for pressure vessels with ellipsoidalheads

6.5 for the pressure vessels with torispherical head. This difference is computed as

Difference(%) =

(CB shell

reference− 1

)× 100 (6.65)

Note that only the collapse loads calculated in [72], which are independent from

elastic properties, are comparable to the CB shell results. For this reason, differences

with [72] for elastic limits and shakedown loads are intentionally omitted in Table

6.5.

Considering the set of models investigated, the following conclusions can be

drawn from the results:

• Regarding the elastic limit (pE)

- The elastic limit decreases as the ratio a/b of the pressure vessels increases;

- The discrepancy with results from literature increases as the wall thickness

increases. The CB shell results are slightly higher (2.2-5.7%) than the results

from [76] for pressure vessels with radius to thickness ratio a/h = 10. The

results are in very good agreement for structures with a/h > 10.

Assuming that all elastic results from solid models [76] are accurate, the rea-

sons for the higher differences observed for models with higher thicknesses can

be:

(a) the limitations of the shell model. When the wall of the pressure vessel

becomes excessively thick, the actual elastic stress state in the knuckle of the

vessel head, where the curvature is higher, is no longer well reproduced by

the shell hypothesis of plane stress state and/or the kinematic assumption of

fibers passing through master nodes remaining straight.

(b) the difference in the way that pressure load is applied. In solid models,

77

Figure 6.10: Examples of collapse mechanisms for pressure vessels with ellipsoidalheads under internal pressure.

78

Figure 6.11: Examples of collapse mechanisms for pressure vessels with torisphericalheads under internal pressure.

Table 6.4: Difference between limit loads calculated with CB shell and results fromYEOM and ROBINSON [76] for pressure vessels with ellipsoidal heads. pE = elasticlimit; pLA = collapse load; pSD = shakedown load

CB Shell vs Reference [76]

Difference(%)a/b a/h pE pLA pSD2 10 — -0.4 —2 25 0.9 0.0 1.72 100 0.1 -0.1 0.42 300 0.9 -0.3 2.1

2.5 10 — -0.6 —2.5 25 0.0 -0.1 5.02.5 100 0.0 -0.1 2.42.5 300 0.0 0.1 3.2

3 10 — -1.0 —3 25 0.9 -0.2 5.23 100 -0.1 0.2 4.43 300 0.1 0.2 3.9

4 10 — -1.5 —4 25 1.2 0.0 4.94 100 -0.1 0.3 4.94 300 0.3 0.4 4.7

the internal pressure is obviously applied in the inner surface of the pressure

vessel. For shell elements the pressure is applied in the shell surface which

coincides with the mean surface of the pressure vessel. For thin-walled pressure

vessels this difference is insignificant. However, when the wall is thick the

total load applied in regions of high curvature in the shell model is higher

79

Table 6.5: Difference between limit loads calculated with CB shell and results fromYEOM and ROBINSON [76] and MAKRODIMOPOULOS [72] for pressure vesselswith torispherical heads. pE = elastic limit; pLA = collapse load; pSD = shakedownload


Difference(%)a/b r/a a/h pE pLA pSD1 1 25 2.0 -0.1 -0.1

1.563 0.2 25 1.5 0.0 0.01.563 0.3 14.29 3.2 -0.1 -0.1

2 0.12 25 1.4 -0.9 4.92 0.2 25 1.1 -0.4 4.82 0.3 14.29 1.8 -0.2 5.12 0.3 25 0.8 -0.3 5.4

2.5 0.12 10 5.7 -1.9 8.52.5 0.12 25 0.9 -1.1 4.92.5 0.12 50 0.2 -0.3 5.32.5 0.2 14.29 1.8 -1.1 5.62.5 0.2 25 0.9 -0.5 4.92.5 0.2 50 0.4 0.0 5.62.5 0.3 10 2.2 -1.8 6.12.5 0.3 14.29 1.4 -0.8 5.42.5 0.3 25 0.6 -0.3 5.22.5 0.3 50 0.1 -0.2 5.8

3.333 0.12 10 4.7 -2.8 7.53.333 0.12 14.29 2.3 -1.6 5.43.333 0.12 25 0.8 -0.6 4.73.333 0.12 50 0.0 0.1 5.73.333 0.2 14.29 1.4 -1.0 4.93.333 0.2 25 0.6 -0.3 4.83.333 0.2 50 0.2 0.3 4.9


Difference(%)R/a r/a a/h pE pLA pSD1.5 0.1 25 — -1.2 —1.5 0.4 25 — 2.0 —

than the total load applied in the respective regions of the solid model. The

effect of this difference can be demonstrated by re-running the CB shell model

corresponding to the torispherical vessel with a/b = 2.5, r/a = 0.12 and

a/h = 10. In this new model we applied a correction factor to the pressure

load to make it compatible to the total load applied in the solid model. As a

result the difference in the elastic limit decreased from 5.7% to 4%.

• Regarding the collapse loads (pLA)

- The instantaneous collapse limit also decreases as the ratio a/b of the pressure

vessel increases;

80

- The CB shell results are very close to the results from references [72, 76].

The maximum discrepancy is only 2.8% in absolute value.

• Regarding the shakedown limit (pSD)

- For the models presenting 2pE < pLA the mechanism of failure is alternate

plasticity and the shakedown load is pSD = 2pE. For the remaining models the

mechanism of failure is incremental collapse and the limit load is pSD ≈ pLA;

- The comparison of shakedown loads shows the larger discrepancies (up to

8.5%) between the CB shell results and the results from [76]. For models

with low a/h, such differences can be partially explained by the same reasons

presented for the elastic limits. Nevertheless, differences can also be originated

by drawbacks faced in [76] when the shakedown limits were calculated by

means of incremental analyses.

Indeed, taking the case of models which mechanism of failure is alternate

plasticity, the shakedown load must be theoretically equal to twice the elastic

limit since the load domain has only two vertices, one of them being zero [84].

Therefore, considering that the computation of the elastic limit is much easier

than calculating the shakedown limit, one can expect that elastic limits from

[76] are more accurate than their shakedown loads. Consequently, differences

observed in shakedown loads should be smaller, i.e. equal to differences obser-

ved for elastic limits. The differences in shakedown loads should be, at least,

closer to the differences in elastic limits assuming, in this case, that elastic

limits from [76] are given in a nodal basis and that some smoothing technique

was used to compute them.

Note that, for some models of ellipsoidal vessels which incremental collapse

is the failure mechanism, the difference between the shakedown limit and the

collapse load calculated in [76] is slightly higher than that observed with the

respective CB shell models. This larger difference may be due to difficulty

encountered in [76] in computing the shakedown load by means of incremental

analysis when the applied load approaches the collapse limit and the stiffness

matrix becomes near singular.

Finally we build an additional model of the torispherical vessel with a/b = 2.5,

r/a = 0.12 and a/h = 10 which presented the maximum differences listed in Ta-

ble 6.5 using the 6-node triangular mixed finite element described in subsection

6.6.1. The objective is to compare the elastic stresses of this solid model with the

corresponding CB shell model.

After a mesh refinement procedure, we obtain the following limit loads from the

solid model: pE = 0.337, pLA = 1.128 and pSD = 0.674. The shakedown mechanism

81

Figure 6.12: Equivalent elastic stress plots from solid (left) and shell (right) modelsfor the pressure vessel with torispherical head. Geometrical parameters of the vessel:a/b = 2.5, r/a = 0.12 and a/h = 10.

of this vessel is alternate plasticity, as correctly predicted by the respective CB shell

model, and the shakedown limit is twice the elastic limit, as expected.

Figure 6.12 shows the plots of von Mises elastic stresses in the knuckle of the

vessel for both solid (left) and shell (right) models. This is the region of higher

stresses in this vessel. The shell model is not taking into account the previously

discussed difference in the way that pressure load is applied.

The maximum von Mises stress computed by the solid model is 7.1% above the

respective stress in shell model. This difference should be 4.6% whether correction in

the applied pressure were considered. Despite this difference, we note from Fig. 6.12

that the shell model can reproduce quite well the stress state in the knuckle with

some slight discontinuities observed when the curvature of the wall changes.

In this case, the difference in the value computed by the solid and the shell model

for the load amplification factor preventing shakedown is acceptable for engineering

purposes.

Note that, as a general rule, care must be taken when using structural elements

whether alternate plasticity is identified as the critical mechanism due to the local

nature of this phenomenon. When there are regions of concentrated stresses in thick-

walled shells, the confirmation of the load amplification factor computed by the shell

model is accomplished by using solid elements associated to a mesh convergence

study. To this end, sub-modeling techniques may be helpful.

82

Chapter 7

Continuum-based

Three-Dimensional Shell Element

In this chapter a six-node triangular shell element is presented for the solution of

3-D shakedown problems. The element is formulated employing the same strategy

used in the development of the beam element and axisymmetric shell elements pre-

sented in the two previous chapters. Hence, the element is built combining the

continuum-based approach [12, 13, 67] and a two-field mixed formulation. In this

case displacements and velocities are quadratic along the element mid-surface and

linear through its thickness. In each layer, stresses are interpolated linearly with

respect to mid-surface parameters and also linearly across the layer thickness. Con-

tinuity of stresses through layers may or may not be enforced, giving two alternative

versions of the finite element.

The proposed CB shell element differs from the family of mixed two-field trian-

gular shell elements of Bathe and co-workers [85–87] where displacement and strain

fields are interpolated and combined in a particular way with the main objective

of avoiding locking, usually present in pure displacement formulations. This shell

element also differs from the mixed (strain-displacement) triangular shell element

originally proposed by Argyris and co-workers [88–91] and further modified by Cor-

radi and Panzeri [62, 63] for application in sequential limit analysis.

In what follows we provide details for implementation of this CB shell element.

Sections are organized in the same way as in the two previous chapters. Note

that some intermediate steps of calculation presented for the CB beam element are

omitted in this chapter to make the text more concise. On the other hand, our

intention is to make the implementation of each finite element easier by giving in its

respective chapters equations in their final form. For this reason and for the sake

of clarity of the exposition, some formulae introduced in previous chapters may be

repeated here.

After the presentation of the element formulation we show the examples of ap-

83

Figure 7.1: Continuum-based shell element: displacement nodes (master and slave)and stress nodes

plication from [92] to assess the accuracy obtained with the use of the proposed 3-D

shell element. The numerical tests include: (i) analyses of thin and thick-walled

straight pipes under combined loads, (ii) the shakedown analysis of a pipe bend

under internal pressure and in-plane bending and (iii) limit analysis of a cylinder-

cylinder intersection subjected to bending and internal pressure. The proposed

element is also used to estimate the external pressure that produces collapse in a

pipe with plain dents with different depths.


The three-dimensional global (fixed) orthonormal reference frame is denoted byR =

ex, ey, ez ≡ ex1 , ex2 , ex3.


The two-dimensional element i has six master nodes and its corresponding directors,

respectively denoted as:

xa = xaex + yaey + zaez, a=1:6 (7.1)

pa = cos θaxex + cos θaye

y + cos θazez, a=1:6. (7.2)

The respective underlying continuum element i has twelve slave nodes:

xα = xαex + yαey + zαez, α=s1:s12. (7.3)

which are also identified by any of the two alternative label sequences (see Fig. 7.1):

s1:s12 ≡ 1−, 1+, 2−, 2+, . . . , 6−, 6+. (7.4)

84

Then, using the above definitions we can relate master and slave nodes by

xa− = xa − ha

2pa, a=1:6 (7.5)

xa+ = xa + ha

2pa, a=1:6 (7.6)

where ha is the shell thickness at node a. Note that (7.5) and (7.6) and the respective

equations for the CB beam element (see (5.5) and (5.6)) are identical, except for the

number of slave nodes (6 and 3 respectively).

We also assume in the 3-D shell element that the fiber through the master node

a and parallel to the director pa moves rigidly. Then, the velocities va, va− and

va+ of master and slave nodes are related by the angular velocity of the director pa

which, in this case, is given by

ωa = ωaxex + ωaye

y + ωazez. (7.7)

Thus, the relation between velocities of master and slave nodes can be recast in any

of the following equivalent forms:

va− = va − ha

2ωa ∧ pa ≡ va + ωa ∧

(xa− − xa

), a=1:6, (7.8)

va+ = va + ha

2ωa ∧ pa ≡ va + ωa ∧

(xa+ − xa

), a=1:6. (7.9)

or in symbolic matrix notation as[va−

va+

]= T a,i

[va

ωa

], (7.10)

where

T a,i :=

[1 3 Λa−

1 3 Λa+

](7.11)

with 1 3 denoting the 3-D identity tensor and Λa+ representing the skew-symmetric

tensor given by Λa+ij =

∑k εijk(x

a+k − xak) (εijk is the alternating symbol). Matrix

Λa− is analogously defined.

Relations (7.10) can be also collected, for a=1:6, in the intrinsic equation (5.12)

vs,i = T ivi (7.12)

by defining the block-diagonal matrix

T i := diag(T 1,i,T 2,i,T 3,i,T 4,i,T 5,i,T 6,i

)(7.13)

85

and

vs,i :=[vs1 vs2 . . . vs11 vs12

]i,T(7.14)

≡[v1− v1+ . . . v6− v6+

]i,T, (7.15)

vi :=[v1 ω1 . . . v6 ω6

]T. (7.16)

Then, (7.12) in global coordinates reads

[vs,i]R

=[T i]R [

vi]R, (7.17)

with

[vs,i]R

=[vs1x vs1y vs1z . . . vs12

x vs12y vs12

z

]R,T(7.18)

=[v1−x v1−

y v1−z . . . v3+

x v3+y v3+

z

]R,T, (7.19)

[vi]R

=[v1x v

1y v

1z ω

1x ω

1y ω

1z . . . v

6x v

6y v

6z ω

6x ω

6y ω

6z

]R,T, (7.20)

and

[T a,i

]R=

1 0 0 0 za− − za ya − ya−

0 1 0 za − za− 0 xa− − xa

0 0 1 ya− − ya xa − xa− 0

1 0 0 0 za+ − za ya − ya+

0 1 0 za − za+ 0 xa+ − xa

0 0 1 ya+ − ya xa − xa+ 0

R

. (7.21)


The parent triangular prism is

Ω := (ξ, η, ζ) ∈ [0, 1]× [0, 1]× [−1, 1], ξ + η ≤ 1, (7.22)

and the geometry mapping in the i-th finite element domain Bi is

Ω 3 ξ = (ξ, η, ζ) 7→ x(ξ)|Bi =∑

α=s1:s12

gα(ξ)xα, (7.23)

86

where the interpolation functions gα(ξ) are the Lagrange interpolation functions for

triangular coordinates. Defining the parameter ψ(ξ, η) := 1− ξ − η, we write:

g1 = 12ψ(2ψ − 1)(1− ζ), g7 = 2ψξ(1− ζ),

g2 = 12ψ(2ψ − 1)(1 + ζ), g8 = 2ψξ(1 + ζ),

g3 = 12ξ(2ξ − 1)(1− ζ), g9 = 2ξη(1− ζ), (7.24)

g4 = 12ξ(2ξ − 1)(1 + ζ), g10 = 2ξη(1 + ζ),

g5 = 12η(2η − 1)(1− ζ), g11 = 2ηψ(1− ζ),

g6 = 12η(2η − 1)(1 + ζ), g12 = 2ηψ(1 + ζ).

The curvilinear coordinate system R(ξ) = ex, ey, ez, used to enforce the plane

stress assumption, is obtained by the method originally proposed by HUGHES [17],

same as that used in [12].

First we define the base vector normal to the lamina as

ez :=x,ξ ∧ x,η‖x,ξ ∧ x,η‖

, (7.25)

where subscripts “, ξ”and “, η”denote the corresponding derivatives.

Then, we define base vectors that lie in the tangent plane

ex :=vaux1 − vaux2

‖vaux1 − vaux2‖ey :=

vaux1 + vaux2

‖vaux1 + vaux2‖, (7.26)

where the auxiliary vectors vaux1 and vaux2 are computed by

vaux1 =x,ξ‖x,ξ‖

+x,η‖x,η‖

vaux2 = ez ∧ vaux1. (7.27)


for vectors and second order tensors (5.32) with

R(ξ) :=

ex · ex ex · ey ex · ez

ey · ex ey · ey ey · ez

ez · ex ez · ey ez · ez

. (7.28)

The Jacobian of the geometry mapping (J(ξ) in (5.34) for α=s1:s12) is given

by x,ξ x,η x,ζ

y,ξ y,η y,ζ

z,ξ z,η z,ζ

R

=∑

α=s1:s12

xαgα,ξ xαgα,η xαgα,ζ

yαgα,ξ yαgα,η yαgα,ζ

zαgα,ξ zαgα,η zαgα,ζ

R

. (7.29)

and the inverse of the Jacobian is given by J−T = [1/ detJ ]CJ , where CJ is the

87

matrix of cofactors of J

CJ =

y,ηz,ζ − y,ζz,η y,ζz,ξ − y,ξz,ζ y,ξz,η − y,ηz,ξx,ζz,η − x,ηz,ζ x,ξz,ζ − x,ζz,ξ x,ηz,ξ − x,ξz,ηx,ηy,ζ − x,ζy,η x,ζy,ξ − x,ξy,ζ x,ξy,η − x,ηy,ξ

R

. (7.30)

Thus, in view of (5.39), the derivatives of the interpolation functions with respect

to spatial coordinates are gα,x

gα,y

gα,z

=1

detJCJ

gα,ξ

gα,η

gα,ζ

(7.31)

for α=s1:s12.


The assumed displacement and velocity fields are then

u(ξ)|Bi =∑

α=s1:s12

gα(ξ)uα v(ξ)|Bi =∑

α=s1:s12

gα(ξ)vα. (7.32)

The equations (7.32) above can also be recast in the compact intrinsic form of (5.43)

u(ξ)|Bi = Nv(ξ)us,i v(ξ)|Bi = Nv(ξ)vs,i,

where, in this case, the vector vs,i, collecting velocities of all slave nodes, is defined

in (7.14) and

us,i =[us1 us2 . . . us11 us12

]i,T(7.33)

≡[u1− u1+ . . . u6− u6+

]i,T, (7.34)

Nv(ξ) :=[g11 3 g21 3 . . . g111 3 g121 3

]. (7.35)

The velocity field written in curvilinear components, i.e. [v(ξ)]R, is given in

terms of the parameters of interpolation written in global components, i.e. [vi]R, by

(5.52)

[v(ξ)]R(ξ) = [Nv(ξ)]RT (ξ)[T i]R [

vi]R.

where the matrices [Nv(ξ)], R(ξ) and[T i]R

of the three-dimensional CB shell

element are analogous to the respective matrices of the elements presented in the

88

previous chapters, i.e., T i is given in (7.13), the linear map Nv which is intrinsic is

[Nv(ξ)] =

g1 0 0 g12 0 0

0 g1 0 . . . 0 g12 0

0 0 g1 0 0 g12

, (7.36)

and the “rotation” matrix is defined as

R(ξ) := diag(R,R,R,R,R,R,R,R,R,R,R,R), (7.37)

with R(ξ) given in (7.28).


We employ the notation for planar tensors introduced in Sect. 5.1.4. Then, assuming

that in the direction perpendicular to the local laminar surface both strain and stress

components are zero we have

d = dxex ⊗ ex + dye

y ⊗ ey + d(xy)

[1√2


)]+

d(yz)

[1√2

(ey ⊗ ez + ez ⊗ ey

)]+ d(xz)

[1√2

(ex ⊗ ez + ez ⊗ ex

)](7.38)

and

σ = σxex ⊗ ex + σye

y ⊗ ey + σ(xy)

[1√2


)]+

σ(yz)

[1√2

(ey ⊗ ez + ez ⊗ ey

)]+ σ(zx)

[1√2

(ex ⊗ ez + ez ⊗ ex

)]. (7.39)

The relevant compatibility equations for the 3-D CB shell element are then

dx = vx,x, (7.40)

dy = vy,y, (7.41)

d(xy) = 1√2

(vx,y + vy,x) , (7.42)

d(yz) = 1√2

(vy,z + vz,y) , (7.43)

d(xz) = 1√2

(vx,z + vz,x) . (7.44)

7.1.5 Computing strain at a generic point

The computation of the infinitesimal strain rate at a generic (Gauss) point of the un-

derlying continuum and in terms of the interpolation parameters [vi]R is performed

89

by using the formula (5.67)



vi]R.

Note that in this equation the vector of strain rates ∇symx v is given in local la-

minar components whereas the vector of nodal velocities [vi]R is written in global

components.

The above equation is calculated by following the same procedure explained for

the CB beam element in Sect. 5.1.5, i.e.,

1. Find the curvilinear coordinate system R(ξk) = ex, ey, ez given by (7.26)

and (7.25) and the rotation matrix R(ξk) using (7.28);

2. Compute, by using

[vs,i]R(ξk)

= RT (ξk)[T i]R [

vi]R

all slave node velocities in the curvilinear system of the generic point;


xα = RT (ξk) xα; (7.45)

4. For α=s1:s12 find (gα,x, gα,y, gα,z);


(xa, ya, za) (from (7.45)) instead of (xa, ya, za). Then, use ∇xgα = J−T∇ξgαwritten for curvilinear coordinates, i.e. analogous to (7.31).

5. Finally, compute the infinitesimal strain tensor

[∇sym

x v(ξk)]R(ξk)

=[∇sym

x Nv(ξk)]R(ξk) [

vs,i]R(ξk)

. (7.46)

Here (7.46) is written, discarding the transversal strain, as

vx,x

vy,y1√2

(vx,y + vy,x)1√2

(vy,z + vz,y)1√2

(vx,z + vz,x)

R(ξk)

=[G1 . . .G12

]R(ξk)

vs1x

vs1y

vs1z...

vs12x

vs12y

vs12z

R(ξk)

(7.47)

90


Gα(ξ) :=

gα,x 0 0

0 gα,x 0

1√2gα,y

1√2gα,x 0

0 1√2gα,z

1√2gα,y

1√2gα,z 0 1√

2gα,x

R

, α=s1:s12. (7.48)

7.2 Stresses

Stresses are interpolated by employing a similar strategy to that used to interpolate

stresses in the beam element (see Sect. 5.2) and in the axisymmetric shell element

(see Sect. 6.2).

The parent domain Ω =⋃`=1:L Ω` is divided to map each element layer. To this

end a sequence of coordinates −1 = ζ0 < ζ1 < · · · < ζL = 1 is selected. Then, each

parent layer is defined as Ω` := (ξ, η, ζ) ∈ (0, 1) × (0, 1) × (ζ`−1, ζ`), ξ + η ≤ 1.The image of these L parent layers through the geometry mapping (7.23) defines

the spatial layers.

The local coordinate ζ ∈ (−1, 1), to be used in the restricted domain of one

single layer `, is defined as

ζ :=ζ − a`b`

, (7.49)

where a` := 12

(ζ` + ζ`−1

)and b` := 1

2

(ζ` − ζ`−1

).

Stress components in each layer ` are interpolated by linear functions in ξ, η and

ζ. Then, each layer has six stress nodes located at corners and

σx(ξ)|Bi =∑j=1:6

tj(ξ)σj|`x , (7.50)

σy(ξ)|Bi =∑j=1:6

tj(ξ)σj|`y , (7.51)

σ(xy)(ξ)|Bi =∑j=1:6

tj(ξ)σj|`(xy), (7.52)

σ(yz)(ξ)|Bi =∑j=1:6

tj(ξ)σj|`(yz), (7.53)

σ(zx)(ξ)|Bi =∑j=1:6

tj(ξ)σj|`(zx), (7.54)

where the interpolation parameters σj|`x , σ

j|`y , σ

j|`(xy), σ

j|`(yz) and σ

j|`(zx) are interpreted as

the stress components of the j-th stress node of layer ` in its own laminar directions.

91

The interpolation functions are

tj(ξ) = 12(ξjξ + ηjη + ψjψ)(1 + ζj ζ) (7.55)

or

tj(ξ) = 12(ξjξ + ηjη + ψjψ)

[1 +

ζj − a`b2`

(ζ − a`)], (7.56)

with (ξj, ηj, ψj) denoting the triangular coordinates of node j. These are the usual

linear functions, just restricted to one layer, where an appropriate transversal co-

ordinate was defined in (7.49). This change of variable is also used to compute

integrals by numerical quadrature.

Finally, the above interpolation of stress is cast a matrix form correspondent to

(5.73) as

[σ(ξ)]R = [Nσ(ξ)][σ`], (7.57)


[Nσ(ξ)] = [t11 5 . . . t61 5] (7.58)


[σ(ξ)]R =[σx σy σ(xy) σ(yz) σ(zx)

]T, (7.59)[

σ`]

=[σ1|` σ2|` σ3|` σ4|` σ5|` σ6|`]T (7.60)


[σt|`]

=[σt|`x σ

t|`y σ

t|`(xy) σ

t|`(yz) σ

t|`(zx)

]T. (7.61)

7.3 Yield function

The yield function f(σ) for the 3-D CB shell element is obtained from the von Mises

yield criterion (5.78) by considering the plane stress condition (σz = 0) as

f(σ) = σ2x + σ2

y − σxσy −3

2

(σ2

(xy) + csσ2(yz) + csσ

2(zx)

)− σ2

Y , (7.62)

The parameter cs denotes a coefficient that can be set equal to 0 or 1 to consider

or not the contribution of the shear stresses to yielding. Therefore, if cs = 0 a

Kirchhoff-Love-type shell theory is simulated. Otherwise, cs = 1 and shear stresses

are considered as in Mindlin-Reissner-type shell elements.

92


∇σf =

2σx − σy2σy − σx

3σ(xy)

3csσ(yz)

3csσ(zx)

and ∇σσf =

2 −1 0 0 0

−1 2 0 0 0

0 0 3 0 0

0 0 0 3cs 0

0 0 0 0 3cs

(7.63)


For obtaining the discrete strain operator for the three-dimensional CB shell element

we firstly define in the discrete setting for each layer ` a work-conjugate strain rate

vector [d`]

:=[d1|` d2|` d3|` d4|` d5|` d6|`

]T, (7.64)

where [dt|`]

=[dt|`x d

t|`y d

t|`(xy) d

t|`(yz) d

t|`(xz)

]T(7.65)


Then, we follow the steps presented for the CB beam element in Sect. 5.4 to

identify the layer discrete strain operator as

B`,i := B`,islave

[T i]R. (7.66)

where

B`,islave :=

∫Ω`

[As1,`,i . . .As12,`,i

]RRT detJ dξ dη dζ

(7.67)

with the auxiliary matrix A defined as


Then, by considering (7.48) and (7.57), the matrix A is computed as

A`,i(ξ) =

t11 5

t21 5

t31 5

t41 5

t51 5

t61 5

[G1 G2 . . . G11 G12

]R. (7.69)

Alternatively,

A`,i(ξ) =[As1,`,i . . .As12,`,i

]R, (7.70)

93

with

Aα,`,i =

t11 5

t21 5

t31 5

t41 5

t51 5

t61 5

[Gα

]R, (7.71)

or, explicitly,

Aα,`,i =

t1gα,x 0 0

0 t1gα,x 0

1√2t1gα,y

1√2t1gα,x 0

0 1√2t1gα,z

1√2t1gα,y

1√2t1gα,z 0 1√

2t1gα,x

......

...

t6gα,x 0 0

0 t6gα,x 0

1√2t6gα,y

1√2t6gα,x 0

0 1√2t6gα,z

1√2t6gα,y

1√2t6gα,z 0 1√

2t6gα,x

R

. (7.72)


of each layer:

Bi =∑`

B`,i. (7.73)



stress fields between layers is enforced. When stress fields are discontinuous between


Bi =

B1,islave...

BL,islave

[T i]R. (7.74)


In this section we compute the discrete form of the elastic operator for the 3-D

CB shell element denoted by E. This matrix is only used in the data preparation


94



(2.20) and (2.32).

Firstly, we define the elastic stress-strain operator E considering the plane stress

constraint σz = 0. Then,

[E]R

=E

1− ν2

1 ν 0 0 0

ν 1 0 0 0

0 0 1− ν 0 0

0 0 0 1− ν 0

0 0 0 0 1− ν

, (7.75)


Hence, by following the demonstration steps presented for the CB beam element

in section 5.4, we conclude that the inverse of the discrete elastic operator for an

element i is computed by assembling the contributions of each layer ` as

(Ei)−1

=∑`

(E`,i)−1

, (7.76)

where(E`,i)−1

denotes the inverse of the elastic relation for the layer ` of the i-th

element defined as

(E`,i)−1

:=

∫Ω`

[Nσ(ξ)]T[E−1

]R[Nσ(ξ)] detJ dξ dη dζ. (7.77)

7.6 Eliminating singularities in rotational stiff-

ness

Since our shell element has six degrees-of-freedom per node we must deal with the

fact that, when the element is flat, its stiffness matrix has singularities associated

with rotations/angular velocities about the directors [12]. This problem is addressed

with a procedure analogous to that employed by BENSON et al. [13]. In this case, a

small rotational stiffness is added to the element stiffness matrix on a node-by-node

basis,

Kaω ←Ka

ω + skpa ⊗ pa, (7.78)

where Kaω is the 3×3 sub-matrix associated with master node a, k is the maximum

value in the sub-matrix diagonal which is multiplied by the small non-dimensional

constant s, and pa is the director of master node a. BENSON et al. [13] suggest a

choice of s within the range [10−6, 10−4]. Here the value s = 10−6 has proven to be

95

an adequate choice in all numerical examples considered.


7.7.1 Straight pipe under combined loading

In the following the 3-D CB shell is used to calculate collapse loads and the shake-

down limits for straight segments of pipes subjected to the combination of bending

moment M and internal pressure p or axial force F .

Limit analysis of a thin-walled pipe

In this first example the objective is to obtain the collapse load diagram for a thin-

walled pipe with closed ends under the combined action of bending moment M and

internal pressure p. The pipe is geometrically defined by its mean radius r and wall

thickness t. The material is elastic-perfectly plastic with yield stress σY .

The analytical solution for this problem considering the von Mises yielding cri-

terion is [93] (M

MC

)2

+

(p

pC

)2

= 1, (7.79)

with the collapse bending moment MC and the collapse pressure pC calculated under

thin-shell hypotheses, i.e.

MC = 4r2tσY , pC =2√3

t

rσY . (7.80)

The finite element approximation to this problem is obtained for a pipe with

σY = 200 MPa and r/t = 25. The length of the model, L = 20r, was determined by

a convergence study to avoid boundary effects.

Due to the geometry of the problem, the pipe is represented by a quarter sym-

metry model. Kinematic boundary conditions are properly imposed to guarantee

the correct structural behavior. Rigid body motion in the direction parallel to the

symmetry planes is avoided by restraining an arbitrary node lying on the symmetry

plane transversal to the pipe.

Bending moment is applied as a linear pressure distribution at the free end cross

section. Internal pressure, p, is applied assuming closed ends condition with an

equivalent axial force equal to pπr2 applied as a constant pressure distribution at

the free end cross section.

The pipe is modeled with 4-layer elements without continuity of stresses between

layers. The mesh, depicted in Fig. 7.2, has 504 elements and 6118 degrees of freedom.

96

Figure 7.2: Finite element mesh for thin-walled pipe (r/t = 25)

Figure 7.3: Collapse load curve of thin-walled pipe (r/t = 25) under internal pressureand bending

The results obtained with the CB shell are compared to the analytical solution

in Fig. 7.3. This comparison demonstrates that the CB shell is able to calculate

limit loads in very close agreement with thin shell theory. For the five points shown

in Fig. 7.3 the maximum relative error is 0.022%. Nevertheless, thin pipes subject

to bending can buckle under smaller loads, mainly due to ovalization of the cross

section.

Limit analysis and shakedown of a thick-walled pipe

In this subsection we analyze a thick walled pipe (r/t = 5) to avoid the large

displacement effects mentioned in the thin-walled pipe example. In this case, the

pipe is geometrically defined by its internal and external radii, ro and ri, respectively.

Firstly we calculate the exact collapse load diagram for a thick-walled pipe under

the combined action of a bending moment M and an axial load F .

Figure 7.4 shows the stress field σx(y) in the pipe cross section when the limit

load is attained considering the hypothesis of elastic-perfectly plastic material. In

97

Figure 7.4: Stress field in the cross section of a thick-walled pipe under maximumcombination of axial tension and bending

this figure we also define the variable y corresponding to the distance between the

axis z and the line separating the part of the cross section in tension from the part

in compression. According to this definition, if y > 0 then F is positive (tension),

and if y < 0 then F is negative (compression).

The collapse axial load is computed as F (y) = 2σY A with A denoting the area

of the section between y = 0 and y = y.

Then, if | y |≤ ri

F (y) = 2σY r2o

[arcsin

(y

ro

)+y√r2o − y2

r2o

]

− 2σY r2i

[arcsin

(y

ri

)+y√r2i − y2

r2i

], (7.81)

otherwise

F (y) = 2σY r2o

[arcsin

(y

ro

)+y√r2o − y2

r2o

]− σY πr2

i . (7.82)

Likewise, for the bending moment, assuming y > 0, we have M (y) = 2σY Q with

Q denoting the first moment of area of the part of the section between y = y and

y = ro.

Thus, if | y |≤ ri then

M (y) = 43σY

[(r2o − y2

)3/2 −(r2i − y2

)3/2], (7.83)

otherwise

M (y) = 43σY(r2o − y2

)3/2. (7.84)

The parametric equations (7.81),(7.82), (7.83) and (7.84) are used to build the

exact collapse curve varying the parameter y between zero and ro. This exact

solution is used to verify the results obtained with the CB shell for a pipe with

98

Figure 7.5: Collapse load curve of thick-walled pipe (r/t = 5) under axial force andbending

ro = 11 mm and ri = 9 mm.

The material properties used in this example are: E = 200 GPa, ν = 0.3,

σY = 200 MPa. The length of the model is the same employed in the previous

example as well as the boundary conditions, the element options and mesh. Bending

moment and axial force are respectively applied at the free end cross section as linear

and constant pressure distributions.

The comparison between exact and finite element results is depicted in Figure

7.5. In this figure axial force and bending moment are normalized by the respective

collapse loads, FC and MC , calculated for thick-walled pipes:

FC = σY π(r2o − r2

i

), MC = 4

3σY(r3o − r3

i

). (7.85)

As can be seen, the CB shell also performs well for this problem. Figure 7.6 shows

the plot of plastic strain rates for this example corresponding to point A in the limit

load diagram (see Fig. 7.5). Figure 7.6 looks like a plot from a full 3-D analysis

because the results are plotted for the underlying continuum elements of the CB

shell.

Finally, we give in Fig. 7.8 the collapse load diagram and the shakedown limit

for the same pipe considering the combined action of internal pressure and bending.

The load domains for both limit and shakedown analyses are shown in Fig. 7.7. Note

in Fig. 7.7 (b) that the load domain for shakedown is the combination of a fixed

(constant) internal pressure and a fully-reversed bending moment.

In the diagrams of Fig. 7.8 the applied bending load is normalized by the limit

99

Figure 7.6: Plastic strain rates for a thick-walled pipe (r/t = 5) under collapse dueto the loads (F/FC ,M/MC) = (0.4924, 0.7162)

Figure 7.7: Load domains for (a) the limit analysis and (b) the shakedown analysisof a thick-walled pipe (r/t = 5) under internal pressure and bending

100

Figure 7.8: Collapse load curve and shakedown limit for a thick-walled pipe (r/t = 5)under internal pressure and bending

load MC given in (7.85). Internal pressure is also normalized by the respective

analytical limit, pC , calculated for a thick-walled straight pipe [50]:

pC =2√3σY ln

(rori

). (7.86)

In the shakedown diagram depicted in Fig. 7.8, cyclic bending loads above seg-

ment AB will produce alternate plasticity whereas cyclic bending loads above seg-

ment BC will cause failure due to incremental collapse.

The load domain in point A consists of the load M varying cyclically between

two extremes of equal magnitude and opposite sign. Therefore, for this point, the

alternate plasticity limit is equal to the elastic limit [84]. The moment ME corres-

ponding to the elastic limit for a thick-walled pipe is

ME =σY π (r4

o − r4i )

4ro. (7.87)

Thus, the exact value of the maximum bending moment preventing alternate plas-

ticity is M exact = 0.71874MC . The limit load obtained with the CB shell is

MFE = 0.71053MC which is 1.2% below the exact solution.

101

Figure 7.9: Pipe bend geometry and loads

7.7.2 Shakedown analysis of pipe bend subjected to internal

pressure and bending moment

Next, we calculate the collapse load diagram and the shakedown limit for a pipe

bend with attached straight pipe sections subjected to the combination of internal

pressure p and in-plane bending moment M (see Fig. 7.9). For shakedown analysis

the applied load consists of constant internal pressure and cyclic opening bending

(0 ≤M ≤Mmax).

This type of structure has been studied by many researchers. A list of references

can be found in the following recent works on this subject: ABDALLA et al. [49],

CHEN et al. [50], TRAN et al. [51]. Most of the previous works on pipe bends are

restricted to limit analysis and just a few focus on shakedown.

The bend is defined by its mean radius r, bend radius of curvature R, and

thickness t. This geometry can be classified by the curvature factor, cf , defined as

cf :=R/r

r/t=Rt

r2. (7.88)

Accordingly, bends with cf ≤ 0.5 are considered highly curved.

The geometrical data and material properties of this example are: r = 10 mm,

t = 2 mm, R = 20 mm, L = 160 mm, E = 200 GPa, ν = 0.3, σY = 200 MPa.

Thus, the present example consists of a thick-walled (r/t = 5) and highly curved

102

Figure 7.10: Finite element mesh using CB shell elements

(cf = 0.4) pipe bend. These properties are taken from an example given in [50].

The finite element model is built, for a quarter of the pipe bend geometry, em-

ploying 6-layer elements without continuity of stresses between layers. Internal pres-

sure assuming closed end condition and bending moment is applied as explained in

the previous example.

The mesh, shown in Fig. 7.10, has 720 elements and 9005 degrees of freedom.

This mesh density is similar to the mesh used in [50], however in the cited work three-

dimensional solid elements were employed with a different technique to perform the

shakedown analysis.

The results obtained with the CB shell are presented in Fig. 7.11. This figure

also shows the collapse load diagram and the shakedown limit found by CHEN et al.

[50]. There are also 5 points from elastic-plastic incremental analyses performed

using Rik´s method in ABAQUS code which were used by these authors to verify

their collapse load results. We remark that all results are normalized by the collapse

limits MC and pC given in (7.85) and (7.86).

It can be seen in Fig. 7.11 that collapse load results obtained with the CB shell

are very close to the results found in [50] and verified with ABAQUS (also in [50]).

The segment BC of the shakedown limit diagram obtained with the CB shell is also

in very good agreement with the results provided in [50]. Cyclic loads above segment

BC cause failure due to incremental collapse. On the other hand, the limit indicated

by segment AB is considerably different from that found in [50]. Cyclic loads above

segment AB produce alternate plasticity.

The determination of the limit against alternate plasticity is much simpler than

the problem of elastic shakedown formulated in Sect. 2.3. Indeed, the formulation

103

Figure 7.11: Collapse load curve and shakedown limit for a thick-walled (r/t = 5)and highly curved (h = 0.4) pipe bend under internal pressure and bending. CL =collapse load; SD = shakedown limit

for alternate plasticity is obtained neglecting the equilibrium constraints so that

yield condition can be checked independently at each point of the body [10, 84, 94].

The accurate determination of the limit preventing alternate plasticity exclusivelly

rely on the accurate computing of elastic stresses.

The load domain represented by point A in the shakedown diagram corresponds

to the bending moment M varying between zero and the maximum load. For this

particular type of load domain the alternate plasticity limit is simply twice the

elastic limit [84]. Then, if the exact elastic solution of this problem were known the

exact limit against alternate plasticity would be easily calculated as the division of

twice the yield stress by the maximum von Mises elastic stress.

When the exact solution is not available, better estimates of the alternate plas-

ticity limit are achieved by improving the finite element elastic results. Then, in

order to check the elastic results produced by the CB shell element, a new model is

built using ANSYS. This model employs quadratic solid elements to represent the

pipe bend and part of the attached straight pipe. The remaining part of the pipe,

far from the region of interest, is modeled with quadratic shell elements. Shell-solid

connection is made using ANSYS multipoint constraint (MPC) feature. The solid

element chosen from ANSYS element library is SOLID186 with displacement for-

mulation. Both full and reduced integration schemes are tested and compared. The

shell element utilized is SHELL 281.

104

Figure 7.12: Shell-solid model for the pipe bend. Detail of the coarser mesh employedin bend region

A mesh solution convergence study was conducted for the part modeled with

solid elements. A coarser mesh was built with four elements through the thickness,

ten around the radius of the bend and twenty around the circumference of the pipe

(see Fig. 7.12). The number of elements in each direction of pipe bend was doubled

and tripled and the maximum von Mises elastic stress for each mesh was computed

to calculate the respective estimates for the load multiplier µ preventing alternate

plasticity. The results of the convergence study are summarized in Table 7.1 which

also provides the result obtained with the CB shell for comparison. It can be seen

that the estimate for the alternate plasticity limit calculated with the CB shell

element is only 2.2% less than the limit computed with the most refined mesh of

ANSYS (mesh 3) considering the full integration scheme or 1.8% less than the limit

calculated with the ANSYS model (mesh 3) for the reduced integration technique.

Figures 7.13 and 7.14 show, respectively, the von Mises stress plot corresponding to

the ANSYS model (full integration) and the CB shell model.

7.7.3 Limit analysis of a cylinder-cylinder intersection

In this example we calculate the collapse load diagram for a cylinder-cylinder inter-

section. The structure is subjected to internal pressure p and a bending moment

M acting on the nozzle. The shell is made of a material with yield stress σY = 234

MPa and is geometrically defined by its internal diameter D = 285 mm and wall

thickness s = 15 mm. The material of the nozzle has yield stress σY = 343 MPa and

its internal diameter and wall thickness are, respectively, d = 20 mm and t = 7.5

105

Figure 7.13: Equivalent elastic stress (von Mises) from ANSYS mesh 3

Figure 7.14: Equivalent elastic stress (von Mises) from CB shell model

106

Table 7.1: Load multiplier µ preventing alternate plasticity for the pipe bend undercyclic bending moment

Maximum von Mises Load multiplierModel Stress [MPa] µ

ANSYS mesh 1 - full integration 771.6 0.5184ANSYS mesh 2 - full integration 799.5 0.5003ANSYS mesh 3 - full integration 796.6 0.5021

ANSYS mesh 1 - reduced integration 801.2 0.4993ANSYS mesh 2 - reduced integration 807.0 0.4956ANSYS mesh 3 - reduced integration 799.8 0.5002

CB shell 814.7 0.4910

Figure 7.15: Cylinder-cylinder intersection: geometry and loads

mm (see Fig. 7.15). The material properties and geometrical data came from the

example provided in [95], however collapse and and shakedown limits for this type

of structure can be also found in the literature for different geometries and loads

[96–100].

Bending moment M is applied using an equivalent pressure distribution at the

free end cross section of the nozzle. Internal pressure, p, is applied assuming closed

end condition, with equivalent axial forces applied as constant pressure distributions

at free end cross sections of the cylinders. Due to the symmetry only a half of the

structure is discretized using the CB shell elements with four layers and without

continuity of stresses between layers. The mesh shown in Fig. 7.16 has 1200 elements

and 14746 degrees of freedom.

A comparison made by VU et al. [95] of the results obtained with their method

with two other results available in literature revealed large differences among all

107

Figure 7.16: Cylinder-cylinder intersection: finite element mesh with CB shell ele-ments

Figure 7.17: Cylinder-cylinder intersection: finite element mesh with ANSYSSHELL93 elements

108

Figure 7.18: Collapse load curve for the cylinder-cylinder intersection

of them. For this reason, in order to have an additional basis for comparison, we

perform some incremental analyses using ANSYS. These models are built employing

quadratic shell elements (SHELL93) from ANSYS library and do not consider large

displacement effects. The mesh used, shown in Fig. 7.17, is defined with the same

number of divisions of the CB shell mesh, resulting in 600 quadrilateral elements.

Finally, Figure 7.18 depicts the comparison of results obtained with the CB shell

and ANSYS. For the sake of completeness, the three results published in [95] are

also plotted. As can be seen, limit loads calculated by direct method using the CB

shell are in very good agreement with incremental analysis results obtained with

ANSYS. Both results are also close to the result 2 provided in [95].

7.7.4 Collapse of Dented Pipelines Subjected to External

Pressure

In this last example we briefly introduce an important problem for offshore pipelines:

the collapse of dented pipes subjected to external pressure. In the numerical tests we

focus on a steel pipe with diameter-to-thickness ratio equal to 18.9. For this value

of D/t the expected failure for a defect-free pipe is driven by plastic collapse rather

than elastic buckling [2, 101]. For dented pipes both phenomena may interact, hence

collapse loads calculated with the proposed element for three different dent depths

are verified against limit loads obtained with non-linear geometric incremental finite

element analyses.

109

In what follows we firstly summarize the significance of dents for the pipeline

integrity. Secondly, we present a short literature review about this subject. After

that we finally carry out the present work.

The objective of this study is to verify whether hypotheses of the classical sha-

kedown theory could be successfully applied to calculate the collapse pressure of a

thick-walled pipe in deep or ultra-deep water depths scenarios. Note, however, that

rather than making an extensive assessment by means of a large set of finite element

models we perform just an initial evaluation considering a small set of models. The

validation of the method for practical purposes depends, of course, on a more com-

prehensive study to be made elsewhere that shall consider different pipe and dent

geometries, material properties, etc.

Significance of dents

A pipeline dent can be formally defined as “a depression which produces a gross

disturbance in the curvature of the pipe wall, caused by contact with a foreign body

resulting in plastic deformation of the pipe wall”[102]. Dents may be classified based

on the geometry of defect. Accordingly, plain dents are those that cause smooth

changes in curvature of the pipe wall, do not contain reduction in wall thickness and

do not change the curvature of adjacent welds. Dents may also be classified based

on the type of constraint. In this case a dent is named as constrained dent if it is

prevented from recover roundness under the influence of internal pressure, otherwise

it is called unconstrained dent.

Dents are a common type of mechanical damage found in pipelines. They repre-

sent a threat to pipeline integrity due to different reasons, namely:

(i) Dents produce stress and strain concentration reducing the fatigue life of the

pipeline.

(ii) Dents locally reduce the pipe diameter, hence they can be barriers to in-line

inspection tools (pigs).

(iii) If associated to loss of pipe wall material, dents can be very dangerous because

this type of defect severely reduces burst resistance of the pipe.

(iv) Dents can locally reduce the collapse strength to external pressure of offshore

pipelines installed in deep water regions and can be the initiators of propaga-

ting buckles [103, 104]apud [101].

The development of methods to assess the integrity of dented pipelines is very

important due to the high cost involved in the repair or substitution of a pipeline

containing a dent. For this reason there is a broad literature covering the structural

110

integrity of pipes containing dents. Studies comprise both numerical and experi-

mental tests.

In 1999 a joint industry project sponsored by sixteen international oil and gas

companies, including PETROBRAS, was set to make a comprehensive literature

review to produce a document specifying the “best” methods for assessing defects

in pipelines. As a result of such a project, in 2003 was released the Pipeline Defect

Assessment Manual (PDAM)[102, 105–107] that contains guidance for the assess-

ment of different types of defect including plain dents, kinked dents, smooth dents

on welds, smooth dents containing gouges and smooth dents containing other types

of defects. The modes of failure covered by PDAM for dented pipes are: burst due

to internal pressure and fatigue.

Literature review on collapse of dented pipes under external pressure

There is no specific method in PDAM for the evaluation of collapse due to external

pressure in dented pipelines but only for defect-free pipes. According to PDAM the

method which presented best agreement with published data was the approach of

DNV-OS-F101[2]. In this method three factors that can lead to collapse are taken

into account: out-of-roundness, circumferential yielding and elastic instability.

The formula for local buckling due to external pressure used in the DNV-OS-

F101 Rules for Pipeline Systems, 2010 edition [2, 108] is is based on the difference

between the external pressure pc and the minimum internal pressure pmin at any

point along the pipeline route, i.e.,

pe − pmin ≤pc

γm · γSC. (7.89)

In equation above γm denotes the material safety factor to account for the uncer-

tainty in material properties and γSC is a safety factor based on the safety class

defined from the consequences of a possible failure.

The collapse capacity pc is determined by

(pc − pel) ·(p2c − p2

p

)= pc · pel · pp · f0 ·

D

t(7.90)

where

pel :=2E

1− ν2

(t

Do

)3

(7.91)

pp :=fy · αfab ·2t

D(7.92)

f0 :=Dmax −Dmin

D(7.93)

111

In the set of equations above pel and pp denotes the elastic and plastic collapse

capacities and f0 represents the pipe out-of-roundness. Equation (7.90) is valid for

f0 ≥ 0.5%.

Other parameters are: αfab denoting the fabrication factor; fy denoting the

yield stress to be used in design; E and ν corresponding respectively to the Young’s

modulus and the Poisson’s ratio; D denoting the nominal outside diameter; Dmax

andDmin corresponding respectively to the greatest and the smallest measured inside

or outside diameters; and t denoting the nominal wall thickness of pipe. The wall

thickness needs to account for the fabrication tolerance and the corrosion allowance,

if relevant.

Regarding the particular evaluation of collapse limits for dented pipes under

external pressure, it can be mentioned the work of KYRIAKIDES and CORONA

[101] who proposed curves to assess the maximum external pressure for pipes con-

taining plain dents deeper than 0.5% of the outside diameter. These curves, which

are based on the work of PARK and KYRIAKIDES [109], are compared in [101] to

some experimental studies complemented by numerical results of denting followed

by collapse. The proposed curves showed good agreement with their experimental

data. Some details of the investigation conducted in [101] are given below.

The experiments considered in [101] consisted of pipes made of stainless steel SS-

304 with nominal D/t ratio varying between 18.9 and 33.6. Specimens with length

about 30D were placed resting in a rigid flat plate, then the dent was produced at

mid-length by means of spherical-shaped indenters. No pressure was applied during

denting phase. An indenter with diameter d = 0.4D was employed to produce most

of the dents. A few tests were performed with an indenter with d = 1.6D. The

maximum dent depth after unloading was less than 0.4D.

After indentation phase the tubes were sealed and placed in a pressure vessel

for the application of external pressure. The pressure was slowly increased until the

observation of a sudden drop in pressure which characterized the limit pressure.

The limit pressures observed in the tests were normalized by the smaller value

between the elastic buckling pressure Pc and the yield pressure Po of each pipe, given

by

Pc =2E

1− ν2

(t

Do

)3

and Po =2σY t

Do

(7.94)

where σY denotes the yield stress. The geometrical parameters are the wall thickness

t and the mean diameter Do of the undeformed pipe. Note that Pc = pel in (7.91)

and Po = pp in (7.92) for αfab = 1.

The normalized pressures were used in [101] to build graphs relating the limit

112

Figure 7.19: Geometric parameters of most deformed cross section of a dent

pressure and the ovalization parameter ∆od of the pipe

∆od :=Dmax −Dmin

Dmax +Dmin

. (7.95)

In the equation above Dmax and Dmin respectively denote the maximum and the

minimum distance across of the convex part of the most deformed cross section of

the pipe, as shown in Fig. 7.19.

The experiments described above were numerically reproduced in [101] by me-

ans of non-linear incremental finite element analyses performed in a commercial FE

code. They also carried out additional simulations employing different diameters

of indenter to assess the influence of the dent shape on the limit pressure. Four

node shell elements with reduced integration and hourglass control were used in

[101] to model the pipes. Contact between the pipe and the indenter as well as the

support was assumed to be frictionless. The pipe was assumed to be made of a

J2-type elastoplastic material with isotropic hardening. The stress-strain response

was approximated using Ramberg-Orgood-type curves. The load steps applied in

the numerical tests corresponded to the loads from experiments. The Rik´s method

was employed in [101] to capture the post critical behavior of the pipes after the

instability developed at maximum external pressure. It was found in [101] a good

correlation between experimental tests and finite element results. Moreover, it was

observed a good correlation between the ovalization parameter and the limit pres-

sure.

KYRIAKIDES and CORONA [101] also compared their experimental limit pres-

sure vs. ovalization curves with the post-critical behavior of pipes with different

degrees of uniform initial ovalization. They concluded that, despite the complex

three-dimensional nature of dents, the most important parameter to the determi-

nation of the collapse pressure of a dented pipe is the ovalization ∆od. Then that

authors proposed design curves for estimating the limit pressure of dented pipes

that correspond to the post-critical response of a finite element model of the pipe

with an uniform initial ovalization of 0.5%.

113

Recent works on collapse limits of dented pipelines under combined loads

Although this is not central to the topic discussed in this example, it is worth to

mention the importance of the determination of collapse limits for different types of

load usually present in offshore (and also onshore) pipelines containing dents or not,

namely: bending moment, axial force and internal pressure. Collapse limits can be

assessed separately for each load or in a combined manner. In particular, the new

pipelines to be installed in ultra-deep water depths may face collapse, caused by

large bending deformation in the sag bend during installation phase [108]. Hence,

in what follows we outline some recent works towards the determination of collapse

limits for other types of load giving a glance on the state of the art, but far from

the purpose of covering all literature on the interaction of multiple loads in dented

pipelines.

• IFLEFEL et al. [110] studied by means of incremental finite element analysis

the capacity of a dented pipe under the combined action of internal pressure

and bending moment. In that work the pipe was modeled with solid elements.

The dent shape was produced with the tube pressurized. The load capacity of

the pipe was investigated for pure bending, pure internal pressure and the com-

bination of bending and pressure loads. Interaction diagrams were provided

for this latter load case.

• BLACHUT e IFLEFEL [111] investigated by means of experiments and finite

element models the collapse caused by bending moment of steel pipes contai-

ning gouged dents. They performed five small scale experiments of machined

pipes with outside-diameter-to-wall-thickness ratio D/t = 40. The dents were

generated with a rigid hemispherical indenter. The dent depths varied from

15% to 23% of the outside diameter of the pipe. The pipes were loaded by ben-

ding moment whilst subjected to constant internal pressure. Bending moment

was increased until the maximum load carrying capacity was achieved. The

load capacity was defined in [111] when well-defined drop of the applied load

was identified in the test. Their numerical results presented a good agreement

in the modeling of denting, however a not so good correlation was found in the

modeling of bending which provided underestimated load vs. rotation curves

in comparison to experimental observations.

• LIMAM et al. [112] employed small scale experiment and incremental finite

element analysis with solid elements to investigate the effect of small plain

dents on bending capacity of pressurized pipes. The experiments were con-

ducted on SS-321 tubes with diameters of 31 mm and D/t = 52. Transverse

dents with depths up to 1.7 times the pipe wall thickness were produced in

114

[112] with a cylindrical-shaped indenter. The pipes were unpressurized during

denting phase, but a paraffin wax plug was employed to produce more loca-

lized dent shapes. During bending phase the paraffin was removed and the

pipes were pressurized to approximately one half of the yield limit pressure.

It was found that the dent tends to significantly reduce bending capacity (by

up to 50%).

• HYUN BAEK et al. [113] used full scale tests and incremental elastic-plastic

finite element models to study the influence of dents on the collapse behavior

of pipes under combined internal pressure and in-plane bending. API 5L X65

pipe with diameter of 762 mm and diameter-to-thickness ratio about 44 were

used to produce dents with depth-to-diameter ratios from 2.5% to 20%. The

dents were produced by hemispherical rods with diameters from 40 mm to 320

mm. A closing or opening in-plane bending load was applied to the dented

pipes pressurized to pressure levels ranging from atmospheric pressure to 16

MPa. Moment-bending angle curves for the dented pipe were obtained from

finite element results for a variety of factors. The load bearing capacity of the

dented pipes under the combined load was evaluated by TES (Twice Elastic

Slope) moments. The opening bending mode showed a higher load bearing

capacity than the closing bending mode under the combined load regardless of

dent depth. The increase in dent depth and internal pressure led to a decrease

of the TES moment regardless of the bending modes.

The present work

We assess the collapse resistance of dented pipes under external pressure by em-

ploying the classical shakedown theory presented in Chapter 2 and the 3-D CB shell

element introduced in this chapter. To the best of our knowledge, this is the first

published work where a finite element procedure based in a direct method is used

with this purpose. This lack of literature might be due to the fact that, depending

on the diameter-to-thickness ratio, pipes under external pressure may become insta-

ble before the plastic capacity of the material is exhausted. Additionally, geometric

changes may play a significant role in the determination of the limit load. Both

situations are not covered by the classical shakedown theory.

We consider a pipe with dents of different depths produced by the action of a

spherical-shaped indenter. Pipe and indenter geometries are based on the thickest

dented pipe investigated in [101].

Pipes with D/t = 18.9 are usually installed in water depths around 2000 m [108].

In such water depths is not common the occurrence of dents produced by the impact

of anchors or trawl gears, such as the cases presented in [114]. However, the trend

115

Figure 7.20: Geometry of the dented pipe model

in oil and gas industry is the development of oil-fields in ultra-deep water depths

(deeper than 1500 m). It is also expected an increase in the use of pipelines to

transport the hydrocarbons directly from the field to shore or into a pipeline grid

[108].

The dent geometries used for the direct assessment of the limit pressures are

imported from incremental non-linear finite element models carried-out with ANSYS

code. We employ incremental analysis not only to simulate the denting process and

to get the dent geometry but also to calculate limit pressures. The limit pressures

obtained from incremental models are used to verify the results computed with the

classical shakedown theory. Both incremental and direct analyses are performed

assuming that the material response is elastic-ideally plastic. Hence, the major

difference between collapse loads calculated in ANSYS to those directly calculated

with the proposed CB shell element are the nonlinear geometric effects which are

not covered in the classical shakedown theory.

Model description

The pipe with external diameter D and wall-thickness t is resting in a flat surface

with length LS. The dent is produced by the action of a spherical-shaped indenter

with diameter DI which is pressed against the pipe wall. Both support and indenter

are hypothetically rigid. Non-dimensional geometric parameters of the model are

given in Table 7.2.

The properties of the elastic-perfectly plastic material are σY = 250 MPa, E =

200 GPa and ν = 0.3.

116

Table 7.2: Non-dimensional geometric parameters of the dented pipe modelD/t DI/D LS/D Lp/D

18.9 0.40 0.99 20

Figure 7.21: Boundary Conditions applied in 3-D models

Incremental analysis

Incremental analyses are performed with ANSYS finite element program version

13.1. Due to symmetry of geometry and loads the pipe is represented by quarter of

the structure. Symmetry boundary conditions (see Fig. 7.21) are properly imposed

to nodes lying on the symmetry planes to guarantee the correct kinematic behavior.

Displacements in axial direction and rotations out of the plane normal to axial

direction are prevented in the pipe end. These boundary conditions are equivalent

to the application of symmetry in the pipe end. A single node of the model is

constrained in vertical direction to prevent free-body movement.

The load steps of the simulations are:

1. The indenter is pressed against the pipe wall until the desired initial depth di.

Three different initial dent depths are considered: di/D = 6.3%, 15.7%, and

25.2%;

2. The indenter is withdrawn allowing the spring-back of the pipe;

3. The external pressure p is gradually increased until the collapse. The applied

pressure load is normalized by the analytical limit pressure

pC = σY ln

(R

R− t

)(7.96)

calculated for a thick-walled straight pipe according to Tresca’s flow theory.

117

1

X

Y

Z

JUN 7 201314:00:58

PLOT NO. 1

ELEMENTS

Figure 7.22: Coarser mesh used in the dented pipe models

The pipe is modeled in ANSYS using SHELL93 elements. The indenter and the

support are modeled with rigid elements TARGE170. Frictionless contact elements

CONTA174 are placed where contact is expected.

Mesh density is defined by a convergence study. Figure 7.22 shows the coarser

mesh used in pipe modeling (msh3d-0) which has 1224 shell elements and Fig. 7.23

depicts the finer mesh (msh3d-2) with 8184 elements. Note that the region of inte-

rest, i.e., where the dent is produced, is covered with denser uniform meshes whereas

the region far from dent is modeled with less elements.

In the incremental analysis the limit pressure may be defined as the applied

pressure when the pressure-deflection curve becomes effectively horizontal, i.e. the

displacement increases rapidly with a very small increase in pressure. Figure 7.24

shows pressure-displacement graphs which confirm the decrease in pipe stiffness

when the applied pressure approaches the collapse limit in the incremental models.

These graphs correspond to the vertical displacement (uy) of the node at the dent

apex during pressurization phase. Models are simulated by employing mesh msh3d-2

and a standard Newton-Raphson algorithm.

Table 7.3 gives the collapse pressures calculated in ANSYS as well as the ma-

ximum dent depth observed at load steps 1 to 3: the initial depth, di, enforced by

the indenter, the dent depth after spring-back, dSB, and the final dent depth, df,

at maximum pressure. We employ the subscript ’inc’ to denote collapse pressures

calculated by means of incremental analysis.

118

1

X

Y

Z

JUN 7 201314:03:20

PLOT NO. 1

ELEMENTS

Figure 7.23: Finer mesh used in the dented pipe models

Figure 7.24: Pressure vs. displacement curve for pipes with R/t = 9.45.

119

Table 7.3: Collapse pressures calculated by incremental analysis.d/D(%)

D/t Initial After SB Final pinc/pC6.3 5.8 6.5 0.864

18.9 15.7 14.8 16.1 0.632

25.2 24.0 25.6 0.525

Direct analysis

Collapse loads are obtained by employing the following procedure:

1. The finite element mesh and boundary conditions for each dented pipe are

obtained by importing the nodal coordinates, connectivities and boundary

conditions of its corresponding ANSYS model;

2. The dent geometry corresponds to the deformed configuration of the ANSYS

model after spring-back;

3. The pipe is modeled with 3-D CB shell elements with 4 layers and no conti-

nuity of stresses between layers. The number of layers was determined by a

convergence study;

4. The maximum amplification factor for the pressure load pC is calculated with

the shakedown algorithm.

The normalized collapse pressures pdir/pC calculated by means of direct analysis

are presented in Table 7.4 for the different diameter-to-thickness ratios D/t and

dent depths d/D measured after spring-back. The column Dif(%) in Table 7.4

corresponds to the difference between incremental results and direct analysis results.

This difference is computed by

Dif(%) =

(pdir

pinc

− 1

)× 100. (7.97)

Note in Table 7.4 that collapse loads calculated with the classical shakedown theory

are in good agreement with results from ANSYS models considering second-order

geometric effects. The collapse loads calculated by the direct method are slightly

higher than the respective loads computed in ANSYS code. The maximum difference

showed less than 6.5% which is acceptable considering the degree of simplification

obtained with the direct assessment of the limit pressure.

Figure 7.25 shows, as an example, the plastic strain rate plot correspondent to

the pipe with a dent depth dSB/D = 24.

120

Table 7.4: Collapse pressures calculated by direct analysis.D/t dSB/R(%) pdir/pC Dif(%)

5.8 0.920 6.4

18.9 14.8 0.668 5.8

24.0 0.551 5.1

Figure 7.25: Example of plastic strain rate distribution for the pipe with a dentdepth dSB/D = 24.

The difference in results between incremental and direct analyses is mostly due

to changes in geometry of the incremental models not captured by the classical

shakedown theory. This is confirmed by a new set of direct analyses carried out

employing the dent geometry corresponding to the last converged step of the incre-

mental models. The results of this new set of models is given in Table 7.5. Note

the decrease in the difference between direct and incremental collapse loads which

confirm the impact of the geometrical changes in the collapse load estimate.

Table 7.5: Collapse pressures calculated by direct analysis considering the finalgeometry from the incremental analysis

D/t df/R(%) pdir/pC Dif(%)

6.5 0.882 2.9

18.9 16.1 0.641 1.4

25.6 0.527 0.5

It is also interesting to compare the collapse loads calculated with incremental

and direct analysis and the design curve obtained with the methodology proposed in

[101]. The design curve is obtained by computing the post-critical behavior of a pipe

with an uniform initial ovalization of 0.5%. The ovalized pipe model is discretized

in ANSYS with the same SHELL93 elements employed in the dented pipe models.

121

Figure 7.26: Maximum pressure vs. ovalization.

The diameter-to-thickness ratio and material properties of the ovalized pipe model

are also taken from the dented pipes models. Figure 7.27 shows the deformed shape

of the ovalized pipe for different stages of pressurization. The unpressurized pipe

is shown in Fig. 7.27(a), the geometry correspondent to the maximum external

pressure is depicted in Fig. 7.27(b) and the deformed shape correspondent to an

external pressure p/pC = 0.17 is given in Fig. 7.27(c). It can not be seen any

remarkable difference between the pipe geometry at zero and maximum pressure.

After the maximum pressure is reached the level of ovalization locally increases and

new equilibrium configurations are found for lower external pressure levels. The

comparison of the collapse loads of the dented pipes and the post-critical behavior

of the ovalized pipe is given in Fig. 7.26. As can be seen the curve from the ovalized

pipe model gives collapse pressures lower than those obtained with the dented pipe

models.

122

(a) p = 0.00pc

ANSYS 13.0PLOT NO. 1

(b) p = 0.89pc


(c) p = 0.17pc


Figure 7.27: Deformed shape of a pipe with an initial ovalization of 0.5% for differentlevels of external pressure.

123

Chapter 8

Conclusion

Continuum-based mixed elements have been devised in this thesis for use in limit

and shakedown analysis of pipes and pressure vessels by means of direct methods

based on the classical shakedown theory.

Direct methods have some relevant advantages as compared to classical incre-

mental plasticity procedures:

• the safety assessment of the structure is made with the sole knowledge of the

load range rather than full load history and this allows to work with more

realistic hypothesis in terms of loading;

• the maximum amplification factor assuring shakedown is calculated without

the need of a theoretically generally infinite number of incremental plasticity

analyses, as a consequence, direct methods are in general less expensive.

The advantages above are of great interest in both academic and industrial points

of view. In fact, it was seen in the literature review (see Sect. 1.2) that the shakedown

theory has been improved by contributions of many researchers with the aim of

relaxing, at least in some degree, the original hypotheses of Melan-Koiter’s theorem

making broader the range of application of the shakedown theory. Moreover, the

development of new numerical tools (or the improvement of the existing ones) for

the effective assessment of complex structures by means of direct methods is still an

objective to be pursued.

This work is focused on the development of new structural finite elements which

are applied only in the framework of the classical shakedown theory. Existing theo-

ries considering for effects not covered by the classical hypotheses are not taken into

account. The theoretical background of the classical shakedown theory was given

in Chapter 2 with the aim of making clear the basic assumptions and limits for the

application of the theory as well as to provide the set of equations which are solved

with the algorithm introduced in Chapter 4.

124

Three different types of finite elements were conceived and developed in the

present work following an increasing degree of complexity, namely: a beam element,

an axisymmetric shell element and a three-dimensional shell element. Details for an

easier implementation of the finite elements were provided in the respective chapter

of each element.

It is worth emphasizing that in many aspects, the requirements for good perfor-

mance of a finite element for shakedown analysis are more demanding than those for

linear elastic analysis or even limit analysis. In fact, the critical mechanism ruling

the computation of shakedown load factors may be:

(i) alternating plasticity, where the peak elastic stress is determinant;

(ii) instantaneous collapse, where the limit analysis scenario is configured, or;

(iii) incremental collapse, where both the ideally elastic stresses (representing loa-

dings) and the residual stresses must be properly approximated.

The proposed structural elements were tested with a set of representative exam-

ples with the aim of confirming that the proposed finite elements can produce ac-

curate solutions for problems involving the features listed above. The examples of

application were described in detail for exploring the different characteristics of each

numerical test. Comparison of results obtained with the proposed elements and the

available analytical solutions or with the results from solid finite element models

showed a very good agreement, which demonstrates the good performance of the

elements.

The finite elements were developed following some guidelines, often referred as

the continuum-based (CB) approach, that have been proven successful in linear

elastic analysis (see e.g. [12]); this may partially justify the good performance

in computing elastic solutions. Moreover, the proposed elements adopt a mixed

(velocity/residual stress) interpolation in each stacked layer – an approach that has

proven effective in shakedown under plane conditions (see [10] and [60]).

In addition, the combination of the aforementioned techniques, which seems to

be novel for collapse and shakedown applications, allowed an easier definition of

plastic admissibility constraints and made straightforward the use of the available

solution algorithm for the shakedown problem. Since the elements are multi-layered,

they can be used with a different material in each layer, allowing, for instance, the

simulation of clad pipes; however, this implementation still needs to be made and

tested.

Based on the results obtained the author believes that the objective of the thesis

was achieved. The proposed structural elements can be a good option for the assess-

ment of collapse and shakedown limits of pipes and pressure vessels in comparison to

125

the use of classical solid elements. Nevertheless, no claim is made here with respect

to adequacy of the proposed elements to other linear or non-linear analyses.

Finally it is worth listing some interesting improvements that can be made as

further extensions of this work:

• The consideration of temperature effects and hardening in the material model;

• The implementation of temperature loads;

• The application of the proposed elements with shakedown theories accounting

for finite displacements.

Note that the third improvement above is actually very important for assessing

collapse and shakedown loads of shell-like structures since those structures are in ge-

neral flexible, thus the classic hypothesis of infinitesimal displacements (and strains)

can often become too restrictive. In this sense, theories such as those proposed by

MAIER and co-workers [115, 116] or by WEICHERT an collaborators [75] could be

improved, implemented and tested with the proposed elements or, alternatively, new

methodologies may be devised for the consideration of large displacement effects.

In this work the problem of large displacements is not addressed as previously

mentioned. The exceptions are (1) the comparison made in the assessment of col-

lapse loads for the thick-walled curved pipe (see Sect. 7.7.2) and (2) the estimate of

collapse loads for a dented pipe (see Sect. 7.7.4). In both examples collapse results

considering large displacements are provided for comparison with limit loads calcula-

ted with the classical shakedown theory framework. It was shown that, under some

conditions, the use of the classical shakedown theory can provide acceptable results

for shell structures. Of course that the definition of such cases rely on the execution

of more detailed tests. However, note that this is a limitation of the classical theory

which has no relation to the accuracy obtained with the proposed elements.

126

Bibliography

[1] ASME Boiler & Pressure Vessel Code, Divison VIII – Rules for Contruction of

Pressure Vessels, Division 2 – Alternative Rules, 2011 ed., 2011.

[2] DNV-OS-F101 Submarine Pipeline Systems, October 2010.

[3] CEN EN 13445 Unfired pressure vessels. Brussels, 2002.

[4] DEBORDES, O., NAYROLES, B. “Sur la theorie et le calcul a l’adaptation

des structures elasto-plastiques”, Journal de Mecanique Theorique et Ap-

pliquee, v. 15, n. 1, pp. 1–53, 1976.

[5] KAMENJARZH, J. Limit Analysis of Solids and Structures. CRC Press, 1996.

[6] KOITER, W. T. “General theorems for elastic-plastic solids”. In: Sneddon,

I. N., Hill, R. (Eds.), Progress in Solids Mechanics, v. 1, North-Holland,

Amsterdam, cap. IV, pp. 165–221, 1960.

[7] NGUYEN, Q. S. Stability and Nonlinear Solid Mechanics. John Wiley & Sons

LTD, 2000.

[8] DE SOUZA NETO, E., PERIC, D., OWEN, D. Computational Methods for

Plasticity: Theory and Applications. Wiley, 2008. ISBN: 9780470694527.

[9] MAIER, G., PASTOR, J., PONTER, A. R. S., et al. “Direct Methods of Limit

and Shakedown Analysis”. In: Milne, I., Ritchie, R. O., Karihaloo, B.

(Eds.), Comprehensive Structural Integrity, Pergamon, pp. 637 – 684, Ox-

ford, 2003. ISBN: 978-0-08-043749-1. doi: DOI:10.1016/B0-08-043749-4/

03059-7.

[10] ZOUAIN, N. “Shakedown and Safety Assessment”. In: Stein, E., de Borst, R.,

Hughes, T. J. R. (Eds.), Encyclopedia of Computational Mechanics, John

Wiley, cap. 9, pp. 291–334, 2004. ISBN: 0-470-84699-2.

[11] AHMAD, S., IRONS, B. M., ZIENKIEWICZ, O. C. “Analysis of thick and

thin shell structures by curved finite elements”, International Journal for

127

Numerical Methods in Engineering, v. 2, n. 3, pp. 419–451, 1970. ISSN:

1097-0207. doi: 10.1002/nme.1620020310.

[12] BELYTSCHKO, T., LIU, W. K., MORAN, B. Nonlinear Finite Elements for

Continua and Structures. New York, Wiley, 2001.

[13] BENSON, D. J., BAZILEVS, Y., HSU, M. C., et al. “Isogeometric shell analy-

sis: The Reissner-Mindlin shell”, Computer Methods in Applied Mechanics

and Engineering, v. 199, pp. 276–289, 2010.

[14] BUECHTER, N., RAMM, E. “Shell theory versus degeneration - a compari-

son in large rotation finite element analysis”, International Journal for


1097-0207. doi: 10.1002/nme.1620340105.

[15] HUGHES, T. J., LIU, W. K. “Nonlinear finite element analysis of shells: Part 1.

three-dimensional shells”, Computer Methods in Applied Mechanics and

Engineering, v. 26, n. 3, pp. 331 – 362, 1981. ISSN: 0045-7825. doi:

10.1016/0045-7825(81)90121-3.

[16] HUGHES, T. J., LIU, W. K. “Nonlinear finite element analysis of shells: Part

2. two-dimensional shells”, Computer Methods in Applied Mechanics and

Engineering, v. 27, n. 2, pp. 167 – 181, 1981. ISSN: 0045-7825. doi:

10.1016/0045-7825(81)90148-1.

[17] HUGHES, T. J. R. The Finite Element Method. Prentice-Hall, Englewood

Cliffs, N. J., 1987.

[18] SIMO, J., FOX, D. “On a stress resultant geometrically exact shell model.

Part I: Formulation and optimal parametrization”, Computer Methods in

Applied Mechanics and Engineering, v. 72, n. 3, pp. 267 – 304, 1989. ISSN:

0045-7825. doi: 10.1016/0045-7825(89)90002-9.

[19] BATHE, K. J. Finite Element Procedures. New Jersey, Prentice Hall, 1996.

[20] HAN, W., REDDY, B. D. Plasticity: Mathematical Theory and Numerical

Analysis, v. 9, Interdisciplinary Applied Mathematics. Springer, 1999.

[21] CHRISTIANSEN, E. “Limit analysis of collapse states”. In: Ciarlet, P. G.,

Lions, J. L. (Eds.), Handbook of Numerical Analysis. Vol. IV - Finite Ele-

ment Methods (Part 2) - Numerical Methods for Solids (Part 2), Elsevier

Science B.V., pp. 193–312, 1996.

128

[22] ZOUAIN, N., BORGES, L. A., SILVEIRA, J. L. “An algorithm for shakedown

analysis with nonlinear yield functions”, Computer Methods in Applied

Mechanics and Engineering, v. 191, n. 23-24, pp. 2463–2481, 2002. ISSN:

0045-7825. DOI: 10.1016/S0045-7825(01)00374-7.

[23] PRAGER, W. “Limit Analysis: The Development of a Concept”, Problems of

Plasticity, pp. 3–24, 1972.

[24] WIECHMANN, K., STEIN, E. “Shape optimization for elasto-plastic defor-

mation under shakedown conditions”, International Journal of Solids and

Structures, v. 43, n. 22-23, pp. 7145 – 7165, 2006. ISSN: 0020-7683. doi:

DOI:10.1016/j.ijsolstr.2006.03.025.

[25] NERY, D. E. S. Estados Limite de Componentes Mecanicos Considerando

Encruamento Cinematico Limitado (in Portuguese). Tese de Doutorado,

Universidade Federal do Rio de Janeiro, 2007.

[26] ABDEL-KARIM, M. “Shakedown of complex structures according to various

hardening rules”, International Journal of Pressure Vessels and Piping,

v. 82, n. 6, pp. 427 – 458, 2005. ISSN: 0308-0161. doi: DOI:10.1016/j.

ijpvp.2005.01.007.

[27] MELAN, E. “Der Spannungszustand eines Mises-Henckychen Kontinuums bei

veraenderlicher Belastung”, Sitzber, Ak. Wiss. Wien IIa, v. 147, pp. 73–

78, 1938.

[28] MELAN, E. “Zur Plastizitat des raumlichen Kontinuums”, Archive of Applied

Mechanics, v. 9, pp. 116–126, 1938. ISSN: 0939-1533. Disponıvel em:

<http://dx.doi.org/10.1007/BF02084409>. 10.1007/BF02084409.

[29] KOITER, W. T. “A new general theorem on shakedown of elastic-plastic struc-

tures”. In: Proceedings of the Koninklijke Nederlandse Akademie van We-

tenschappen B.59, pp. 24–34, 1956.

[30] KONIG, J. A., MAIER, G. “Shakedown analysis of elastoplastic structures:

A review of recent developments”, Nuclear Engineering and Design, pp.

81–95, 1981.

[31] WEICHERT, D., MAIER, G. Inelastic Analysis of Structures under Variable

Loads. Kluwer Academic Publishers, 2000.

[32] MROZ, Z., WEICHERT, D., DOROSZ, S. Inelastic Behaviour of Structures

under Variable Loads. Kluwer Academic Publishers, 1995.

129

http://dx.doi.org/10.1007/BF02084409

[33] PYCKO, S., MAIER, G. “Shakedown theorems for some classes of nonassoci-

ative hardening elastic-plastic material models”, International Journal of

Plasticity, v. 11, n. 4, pp. 367–395, 1995.

[34] NGUYEN, Q. S. “On Shakedown Analysis in Hardening Plasticity”, Journal

of the Mechanics and Physics of Solids, 2003.

[35] STAAT, M., (EDS.), M. H., HEITZER, M., et al. “Numerical Methods for Li-

mit and Shakedown Analysis - Deterministic and Probabilistic Problems”.

2003.

[36] STAAT, M., HEITZER, M. “LISA – a European project for FEM-based limit

and shakedown analysis”, Nuclear Engineering and Design, v. 206, n. 2-3,

pp. 151 – 166, 2001. ISSN: 0029-5493. doi: DOI:10.1016/S0029-5493(00)

00415-5.

[37] LANG, H., WIRTZ, K., HEITZER, M., et al. “Cyclic plastic deformation

tests to verify FEM-based shakedown analyses”, Nuclear Engineering and

Design, v. 206, n. 2-3, pp. 235 – 247, 2001. ISSN: 0029-5493. doi: DOI:

10.1016/S0029-5493(00)00438-6.

[38] HEITZER, M. “Plastic limit loads of defective pipes under combined internal

pressure and axial tension”, International Journal of Mechanical Sciences,

v. 44, n. 6, pp. 1219 – 1224, 2002. ISSN: 0020-7403. doi: DOI:10.1016/

S0020-7403(02)00025-5.

[39] VU, D. K., YAN, A. M., NGUYEN-DANG, H. “A primal-dual algorithm for

shakedown analysis of structures”, Computer Methods in Applied Mecha-

nics and Engineering, v. 193, n. 42-44, pp. 4663 – 4674, 2004. ISSN:

0045-7825. doi: DOI:10.1016/j.cma.2004.03.011.

[40] KHOI, V. D., YAN, A.-M., HUNG, N.-D. “A dual form for discretized kinematic

formulation in shakedown analysis”, International Journal of Solids and

Structures, v. 41, n. 1, pp. 267 – 277, 2004. ISSN: 0020-7683. doi: DOI:

10.1016/j.ijsolstr.2003.08.013.

[41] STAAT, M., HEITZER, M., LANG, H., et al. “Direct finite element route

for design-by-analysis of pressure components”, International Journal of

Pressure Vessels and Piping, v. 82, n. 1, pp. 61 – 67, 2005. ISSN: 0308-

0161. doi: 10.1016/j.ijpvp.2004.04.006.

[42] LIU, Y. H., ZHANG, X. F., CEN, Z. Z. “Lower bound shakedown analysis by

the symmetric Galerkin boundary element method”, International Jour-

130

nal of Plasticity, v. 21, n. 1, pp. 21 – 42, 2005. ISSN: 0749-6419. doi:

DOI:10.1016/j.ijplas.2004.01.003.

[43] CHEN, S., LIU, Y., CEN, Z. “Lower bound shakedown analysis by using the

element free Galerkin method and non-linear programming”, Computer

Methods in Applied Mechanics and Engineering, v. 197, n. 45-48, pp. 3911

– 3921, 2008. ISSN: 0045-7825. doi: DOI:10.1016/j.cma.2008.03.009.

[44] BELYTSCHKO, T. “Plane Stress Shakedown Analysis by Finite Elements”,

International Journal of Mechanical Sciences, v. 14, pp. 619–625, 1972.

[45] ABDALLA, H. F., MEGAHED, M. M., YOUNAN, M. Y. “A simplified tech-

nique for shakedown limit load determination of a large square plate with

a small central hole under cyclic biaxial loading”, Nuclear Engineering

and Design, v. 241, n. 3, pp. 657 – 665, 2011. ISSN: 0029-5493. doi:

DOI:10.1016/j.nucengdes.2010.11.019. The International Conference on

Structural Mechanics in Reactor Technology (SMiRT19) Special Section.

[46] ABOU-HANNA, J., MCGREEVY, T. E. “A simplified ratcheting limit method

based on limit analysis using modified yield surface”, International Jour-

nal of Pressure Vessels and Piping, v. 88, n. 1, pp. 11 – 18, 2011. ISSN:

0308-0161. doi: DOI:10.1016/j.ijpvp.2010.12.001.

[47] NGO, N., TIN-LOI, F. “Shakedown analysis using the p-adaptive finite element

method and linear programming”, Engineering Structures, v. 29, n. 1,

pp. 46 – 56, 2007. ISSN: 0141-0296. doi: DOI:10.1016/j.engstruct.2006.

03.033.

[48] ABDALLA, H. F., MEGAHED, M. M., YOUNAN, M. Y. “A simplified tech-

nique for shakedown limit load determination”, Nuclear Engineering and

Design, v. 237, n. 12-13, pp. 1231 – 1240, 2007. ISSN: 0029-5493. doi:

DOI:10.1016/j.nucengdes.2006.09.033. 18th International Conference on

Structural Mechanics in Nuclear Engineering, 18th International Confe-

rence on Structural Mechanics in Nuclear Engineering.

[49] ABDALLA, H. F., MEGAHED, M. M., YOUNAN, M. Y. A. “Shakedown

Limit Loads for 90 Degree Scheduled Pipe Bends Subjected to Steady In-

ternal Pressure and Cyclic Bending Moments”, Journal of Pressure Vessel

Technology, v. 133, n. 3, pp. 031207, 2011. doi: 10.1115/1.4002055.

[50] CHEN, H., URE, J., LI, T., et al. “Shakedown and limit analysis of 90 pipe

bends under internal pressure, cyclic in-plane bending and cyclic thermal

131

loading”, International Journal of Pressure Vessels and Piping, v. 88,

pp. 213 – 222, 2011. ISSN: 0308-0161. doi: 10.1016/j.ijpvp.2011.05.003.

[51] TRAN, T. N., KREISSIG, R., VU, D. K., et al. “Upper bound limit and shake-

down analysis of shells using the exact Ilyushin yield surface”, Computers

& Structures, v. 86, n. 17-18, pp. 1683 – 1695, 2008. ISSN: 0045-7949.

doi: DOI:10.1016/j.compstruc.2008.02.005.

[52] CHEN, H. F., SHU, D. “Lower and upper bound limit analyses for pipeline with

multi-slots of various configurations”, International Journal of Pressure

Vessels and Piping, v. 77, n. 1, pp. 17 – 25, 2000. ISSN: 0308-0161. doi:

DOI:10.1016/S0308-0161(99)00074-5.

[53] CHEN, H. F., PONTER, A. R. S. “Shakedown and limit analyses for 3-D

structures using the linear matching method”, International Journal of

Pressure Vessels and Piping, v. 78, n. 6, pp. 443 – 451, 2001. ISSN:

0308-0161. doi: DOI:10.1016/S0308-0161(01)00052-7.

[54] CHEN, H. F., SHU, D. “Simplified limit analysis of pipelines with multi-

defects”, Engineering Structures, v. 23, n. 2, pp. 207 – 213, 2001. ISSN:

0141-0296. doi: DOI:10.1016/S0141-0296(00)00009-2.

[55] STAAT, M., VU, D. K. “Limit loads of circumferentially flawed pipes and cylin-

drical vessels under internal pressure”, International Journal of Pressure

Vessels and Piping, v. 83, n. 3, pp. 188 – 196, 2006. ISSN: 0308-0161. doi:

DOI:10.1016/j.ijpvp.2006.01.002.

[56] STAAT, M., VU, D. K. “Limit analysis of flaws in pressurized pipes and

cylindrical vessels. Part I: Axial defects”, Engineering Fracture Mechanics,

v. 74, n. 3, pp. 431 – 450, 2007. ISSN: 0013-7944. doi: DOI:10.1016/j.

engfracmech.2006.04.031.

[57] HODGE JR., P. G. Limit Analysis of Rotationally Symmetric Plates and Shells.

Prentice Hall, 1963.

[58] GILL, S. A. The Stress Analysis of Pressure Vessels and Pressure Vessels

Components. Pergamon Press, 1970.

[59] ZOUAIN, N. Analise limite de cascas via otimizacao ( in Portuguese). Tese de

Doutorado, Universidade Federal do Rio de Janeiro, 1982.

[60] ZOUAIN, N., HERSKOVITS, J., BORGES, L. A., et al. “An Iterative Algo-

rithm for Limit Analysis with Non-Linear Yield Functions”, ijss, v. 30,

132

n. 10, pp. 1397–1417, 1993. ISSN: 0020-7683. DOI: 10.1016/0020-

7683(93)90220-2.

[61] ZIENKIEWICZ, O. C., TAYLOR, R. L., ZHU, J. Z. The Finite Element

Method: Its Basis and Fundamentals. Elsevier, 2005.

[62] CORRADI, L., PANZERI, N. “Post-collapse analysis of plates and shells based

on a rigid-plastic version of the TRIC element”, Computer Methods in

Applied Mechanics and Engineering, v. 192, n. 33 - 34, pp. 3747 – 3775,

2003. ISSN: 0045-7825. doi: 10.1016/S0045-7825(03)00373-6.

[63] CORRADI, L., PANZERI, N. “A triangular finite element for sequential limit

analysis of shells”, Advances in Engineering Software, v. 35, n. 10 - 11,

pp. 633 – 643, 2004. ISSN: 0965-9978. doi: 10.1016/j.advengsoft.2004.03.

014. Engineering Computational Technology.

[64] ROCKAFELLAR, R. T. Convex Analysis. Princeton University Press, 1970.

[65] HALPHEN, B., NGUYEN, Q. S. “Sur les materiaux standard generalises”, J.

Mecanique, pp. 1–37, 1975.

[66] BORGES, L. A., ZOUAIN, N., HUESPE, A. E. “A Nonlinear Optimization

Procedure for Limit Analysis”, ejma, v. 15, n. 3, pp. 487–512, 1996.

[67] CURIEL SOSA, J. L., GIL, A. J. “Analysis of a continuum-based beam element

in the framework of explicit-FEM”, Finite Elements in Analysis and De-

sign, v. 45, pp. 583–591, 2009.

[68] MARTINS, R. R., ZOUAIN, N., BORGES, L. “A mixed formulation for a

continuum-based beam element”. In: COBEM 21st International Con-

gress of Mechanical Engineering, Natal, RN, Brazil, 2011.

[69] KONIG, J. A. Shakedown of elastic-plastic structures. Elsevier, Amsterdam,

1987.

[70] FRANCO, J. R. Q., PONTER, A. R. S. “A general approximate technique

for the finite element shakedown and limit analysis of axisymmetrical

shells. Part 2: Numerical applications”, International Journal for Nume-

rical Methods in Engineering, v. 40, n. 19, pp. 3515–3536, 1997. ISSN:

1097-0207.

[71] HACHEMI, A., WEICHERT, D. “Numerical shakedown analysis of damaged

structures”, Computer Methods in Applied Mechanics and Engineering,

v. 160, pp. 57 – 70, 1998. ISSN: 0045-7825. doi: 10.1016/S0045-7825(97)

00283-1.

133

[72] MAKRODIMOPOULOS, A. “Computational formulation of shakedown analy-

sis as a conic quadratic optimization problem”, Mechanics Research Com-

munications, v. 33, n. 1, pp. 72 – 83, 2006. ISSN: 0093-6413. doi:

10.1016/j.mechrescom.2004.10.014.

[73] MUSCAT, M., MACKENZIE, D., HAMILTON, R. “Evaluating shakedown

under proportional loading by non-linear static analysis”, Computers &

Structures, v. 81, n. 17, pp. 1727 – 1737, 2003. ISSN: 0045-7949. doi:

10.1016/S0045-7949(03)00181-0.

[74] VU, D. K., STAAT, M. “Shakedown analysis of structures made of materials

with temperature-dependent yield stress”, International Journal of Solids

and Structures, v. 44, n. 13, pp. 4524 – 4540, 2007. ISSN: 0020-7683. doi:

10.1016/j.ijsolstr.2006.11.038.

[75] WEICHERT, D., HACHEMI, A. “Influence of geometrical nonlinearities on

the shakedown of damaged structures”, International Journal of Plas-

ticity, v. 14, n. 9, pp. 891 – 907, 1998. ISSN: 0749-6419. doi:

10.1016/S0749-6419(98)00035-7.

[76] YEOM, D., ROBINSON, M. “Numerical analysis of the elastic-plastic behavi-

our of pressure vessels with ellipsoidal and torispherical heads”, Interna-

tional Journal of Pressure Vessels and Piping, v. 65, n. 2, pp. 147 – 156,

1996. ISSN: 0308-0161. doi: 10.1016/0308-0161(94)00174-H.

[77] MARTINS, R. R., ZOUAIN, N., BORGES, L. A., et al. “A Continuum-Based

Mixed Axisymmetric Shell Element for Limit and Shakedown Analysis”,

Journal of the Brazilian Society of Mechanical Sciences and Engineering,

pp. 1–20, 2013. ISSN: 1678-5878. doi: 10.1007/s40430-013-0059-1. Dis-

ponıvel em: <http://dx.doi.org/10.1007/s40430-013-0059-1>.

[78] DEMIR, H. H. “Cylindrical shells under ring loads”, Journal of the Structural

Division, Proc. of the A.S.C.E., v. 91, n. 3, pp. 71 – 98, 1965.

[79] SAWCZUK, A., HODGE, JR., P. G. “Comparison of yield conditions for cir-

cular cylindrical shells”, Journal of the Franklin Institute, v. 269, n. 5,

pp. 362 – 374, 1960. ISSN: 0016-0032. doi: 10.1016/0016-0032(60)90172-1.

[80] CHAKRABARTY, J. Applied Plasticity, Second Edition. Springer, 2010. ISBN:

978-0-387-77674-3. doi: 10.1007/978-0-387-77674-3.

[81] FRANCO, J. R. Q., BARROS, F. B., PONTER, A. R. “An improved adaptive

formulation for the computation of limit analysis problems of thin axisym-

134

http://dx.doi.org/10.1007/s40430-013-0059-1

metrical shells”, Computers & Structures, v. 87, n. 23 - 24, pp. 1602 – 1610,

2009. ISSN: 0045-7949. doi: 10.1016/j.compstruc.2009.08.014.

[82] KRABBENHØFT, K., LYAMIN, A. V., SLOAN, S. W., et al. “An interior-

point algorithm for elastoplasticity”, International Journal for Numerical

Methods in Engineering, v. 69, n. 3, pp. 592–626, 2007. ISSN: 1097-0207.

doi: 10.1002/nme.1771.

[83] TIMOSHENKO, S., WOINOWSKY-KRIEGER, S. Theory of Plates and Shells.

McGraw-Hill, 1959.

[84] KRABBENHØFT, K., LYAMIN, A., SLOAN, S. “Bounds to Shakedown Loads

for a Class of Deviatoric Plasticity Models”, Computational Mechanics,

v. 39, pp. 879–888, 2007. ISSN: 0178-7675. 10.1007/s00466-006-0076-3.

[85] BATHE, K.-J., LEE, P.-S. “Measuring the convergence behavior of shell analy-

sis schemes”, Computers & Structures, v. 89, n. 3-4, pp. 285 – 301, 2011.

ISSN: 0045-7949. doi: 10.1016/j.compstruc.2010.09.007.

[86] KIM, D.-N., BATHE, K.-J. “A triangular six-node shell element”, Computers

& Structures, v. 87, n. 23 - 24, pp. 1451 – 1460, 2009. ISSN: 0045-7949.

doi: 10.1016/j.compstruc.2009.05.002.

[87] LEE, P.-S., BATHE, K.-J. “Development of MITC isotropic triangular shell

finite elements”, Computers & Structures, v. 82, n. 11-12, pp. 945 – 962,

2004. ISSN: 0045-7949. doi: 10.1016/j.compstruc.2004.02.004.

[88] ARGYRIS, J., PAPADRAKAKIS, M., APOSTOLOPOULOU, C., et al. “The

TRIC shell element: theoretical and numerical investigation”, Computer

Methods in Applied Mechanics and Engineering, v. 182, n. 1 - 2, pp. 217

– 245, 2000. ISSN: 0045-7825. doi: 10.1016/S0045-7825(99)00094-8.

[89] ARGYRIS, J., PAPADRAKAKIS, M., KARAPITTA, L. “Elasto-plastic analy-

sis of shells with the triangular element TRIC”, Computer Methods in

Applied Mechanics and Engineering, v. 191, n. 33, pp. 3613 – 3636, 2002.

ISSN: 0045-7825. doi: 10.1016/S0045-7825(02)00308-0.

[90] ARGYRIS, J., PAPADRAKAKIS, M., MOUROUTIS, Z. S. “Nonlinear dy-

namic analysis of shells with the triangular element TRIC”, Computer

Methods in Applied Mechanics and Engineering, v. 192, n. 26 - 27, pp. 3005

– 3038, 2003. ISSN: 0045-7825. doi: 10.1016/S0045-7825(03)00315-3.

[91] ARGYRIS, J., TENEK, L., PAPADRAKAKIS, M., et al. “Postbuckling per-

formance of the TRIC natural mode triangular element for Isotropic and

135

laminated composite shells”, Computer Methods in Applied Mechanics

and Engineering, v. 166, n. 3 - 4, pp. 211 – 231, 1998. ISSN: 0045-7825.

doi: 10.1016/S0045-7825(98)00071-1.

[92] MARTINS, R. R., ZOUAIN, N., BORGES, L. A., et al. “A Continuum-Based

Mixed Shell Element for Shakedown Analysis”, p. Submitted, 2013.

[93] YAN, A. M., NGUYEN-DANG, H. “Direct finite element kinematical approa-

ches in limit and shakedown analysis of shells and elbows”. In: Weichert,

D., Maier, G. (Eds.), Inelastic Analysis of Structures under Variable Lo-

ads, pp. 233–254. Kluwer Academic Publishers, Dordrecht, 2000.

[94] POLIZZOTO, C. “On the conditions to prevent plastic shakedown of structures:

Part II - The plastic shakedown limit load”, Journal of Applied Mechanics,

, n. 60, pp. 20–25, 1993.

[95] VU, D. K., STAAT, M., TRAN, I. T. “Analysis of pressure equipment by

application of the primal-dual theory of shakedown”, Communications in


1099-0887. doi: 10.1002/cnm.891.

[96] SIMON, J.-W., WEICHERT, D. “Numerical lower bound shakedown analysis

of engineering structures”, Computer Methods in Applied Mechanics and

Engineering, v. 200, n. 41-44, pp. 2828 – 2839, 2011. ISSN: 0045-7825.

doi: 10.1016/j.cma.2011.05.006.

[97] STAAT, M., HEITZER, M. “LISA - a European project for FEM-based limit

and shakedown analysis”, Nuclear Engineering and Design, v. 206, n. 2

- 3, pp. 151 – 166, 2001. ISSN: 0029-5493. doi: 10.1016/S0029-5493(00)

00415-5.

[98] STAAT, M., HEITZER, M. Numerical methods for limit and shakedown analy-

sis: deterministic and probabilistic problems. Relatorio tecnico, John von

Neumann Institute for Computing, 2003.

[99] TRAN, T. N., KREISSIG, R., STAAT, M. “Probabilistic limit and shakedown

analysis of thin plates and shells”, Structural Safety, v. 31, n. 1, pp. 1 –

18, 2009. ISSN: 0167-4730. doi: 10.1016/j.strusafe.2007.10.003.

[100] VU, D., YAN, A., NGUYEN-DANG, H. “A primal-dual algorithm for sha-

kedown analysis of structures”, Computer Methods in Applied Mechanics

and Engineering, v. 193, n. 42–44, pp. 4663 – 4674, 2004. ISSN: 0045-7825.

doi: 10.1016/j.cma.2004.03.011.

136

[101] KYRIAKIDES, S., CORONA, E. “6 - Collapse of Dented Pipes Under

External Pressure”. In: Mechanics of Offshore Pipelines, Elsevier Sci-

ence Ltd, pp. 164 – 180, Oxford, 2007. ISBN: 978-0-08-046732-0. doi:

10.1016/B978-008046732-0/50006-4.

[102] COSHAM, A., HOPKINS, P. “The effect of dents in pipelines – guidance in

the pipeline defect assessment manual”, International Journal of Pressure

Vessels and Piping, v. 81, n. 2, pp. 127 – 139, 2004. ISSN: 0308-0161.

doi: 10.1016/j.ijpvp.2003.11.004. The 10th International Conference on

Pressure Vessel Technology (ICPVT-10).

[103] DEMARS, K. R. “Pipeline failure: A need for improved analyses and site

surveys”. In: Proc. Offshore Technology Conference 4, n. OTC2966, pp.

63–70, 1977.

[104] STRATING, J. “A survey of pipelines in the North Sea. Incidents during

installation, testing and operation”. In: Proc. Offshore Technology Con-

ference 3, n. OTC4069, pp. 25–32, 1981.

[105] COSHAM, A., HOPKINS, P. The Pipeline Defect Assessment Manual - A

Report to the PDAM Joint Industry Project - Confidential, May 2003.

[106] MACDONALD, K., COSHAM, A. “Best practice for the assessment of defects

in pipelines – gouges and dents”, Engineering Failure Analysis, v. 12, n. 5,

pp. 720 – 745, 2005. ISSN: 1350-6307. doi: 10.1016/j.engfailanal.2004.12.

011.

[107] COSHAM, A., HOPKINS, P., MACDONALD, K. “Best practice for the

assessment of defects in pipelines – Corrosion”, Engineering Failure

Analysis, v. 14, n. 7, pp. 1245 – 1265, 2007. ISSN: 1350-6307. doi:

10.1016/j.engfailanal.2006.11.035.

[108] OLAV FYRILEIV, OLAV AAMLID, A. V., COLLBERG, L. “Deepwater

pipelines – status, challenges and future trends”. In: Proceedings of the

Institution of Mechanical Engineers, Part M:, July 2012. doi: 10.1177/

1475090212450659.

[109] PARK, T., KYRIAKIDES, S. “On the collapse of dented cylinders under

external pressure”, International Journal of Mechanical Sciences, v. 38,

n. 5, pp. 557 – 578, 1996. ISSN: 0020-7403. doi: 10.1016/0020-7403(95)

00065-8.

137

[110] IFLEFEL, I., MOFFAT, D., MISTRY, J. “The interaction of pressure and

bending on a dented pipe”, International Journal of Pressure Vessels and

Piping, v. 82, n. 10, pp. 761 – 769, 2005. ISSN: 0308-0161. doi: 10.1016/

j.ijpvp.2005.06.002.

[111] BLACHUT, J., IFLEFEL, I. “Collapse of pipes with plain or gouged dents by

bending moment”, International Journal of Pressure Vessels and Piping,

v. 84, n. 9, pp. 560 – 571, 2007. ISSN: 0308-0161. doi: 10.1016/j.ijpvp.

2007.04.007.

[112] LIMAM, A., LEE, L.-H., KYRIAKIDES, S. “On the collapse of dented tubes

under combined bending and internal pressure”, International Journal of

Mechanical Sciences, v. 55, n. 1, pp. 1 – 12, 2012. ISSN: 0020-7403. doi:

10.1016/j.ijmecsci.2011.10.005.

[113] HYUN BAEK, J., PYO KIM, Y., SIK KIM, W., et al. “Load bearing capacity

of API X65 pipe with dent defect under internal pressure and in-plane

bending”, Materials Science and Engineering: A, v. 540, n. 0, pp. 70 –

82, 2012. ISSN: 0921-5093. doi: 10.1016/j.msea.2012.01.078.

[114] MACDONALD, K., COSHAM, A., ALEXANDER, C., et al. “Assessing me-

chanical damage in offshore pipelines – Two case studies”, Engineering

Failure Analysis, v. 14, n. 8, pp. 1667 – 1679, 2007. ISSN: 1350-6307. doi:

10.1016/j.engfailanal.2006.11.074. Papers presented at the Second Inter-

national Conference on Engineering Failure Analysis (Toronto, Canada,

12-15 September 2006) Part II.

[115] MAIER, G., GRIERSON, D., BEST, M. “Mathematical programming

methods for deformation analysis at plastic collapse”, Computers &

Structures, v. 7, n. 5, pp. 599 – 612, 1977. ISSN: 0045-7949. doi:

http://dx.doi.org/10.1016/0045-7949(77)90002-5.

[116] ARDITO, R., COCCHETTI, G., MAIER, G. “On structural safety assess-

ment by load factor maximization in piecewise linear plasticity”, European

Journal of Mechanics - A/Solids, v. 27, n. 5, pp. 859 – 881, 2008. ISSN:

0997-7538. doi: http://dx.doi.org/10.1016/j.euromechsol.2007.11.011.

138

Documents

CONTINUUM-BASED FINITE ELEMENTS FOR SHAKEDOWN …w2.files.scire.net.br/atrio/ufrj-pem_upl/THESIS/... · the computation of collapse loads and shakedown limits of pipelines and pressure