w.
. ...fi...

1..

.. r ...
o. .. i. a
tn... «1.. z.
. .. t. .
.74.!) .4? “....

.Jx

fl.” ...m.

.é...._..§._

 

 

 

 

 

.1 ...
(u

. 1.5..

.151:
:

 

 

I
an
5. .....s:
...... z .

 

 

 

,.. .. . .: A , ... .., 27.1.3111}:H . . ‘ . _.::.:c. I; , _. .. . .. ,,. .. ..u .
n.,.....~.a...w._n.,._% ...?u.~...u.a.§e.m§§ 4%.? : . 4.... imfiqfifig . . ...w. 5%{3 ax...

, x I . I. :
.. ... .3. ... ... . . . :4...

 

LIBRARY
A Michigan State
C“ -. University
"Q Q} 3“)
(«f/116 3 a

This is to certify that the
dissertation entitled

PSEUDOELASTIC DEFORMATION OF SHAPE MEMORY
ANNULAR PLATES SUBJECT TO NORMAL TRACTIONS

presented by

Yuwei Chi

has been accepted towards fulfillment
of the requirements for the

Ph. D. degree in Engineering Mechanics

 

 

, %'

 

Major Wssor’s Signature

4—” 2’ 05/

Date

MSU is an Affirmative Action/Equal Opportunity Institution

 

PLACE IN RETURN Box to remove this checkout from your record.
To AVOID FINES return on or before date due.
MAY BE RECALLED with earlier due date if requested.

 

DATE DUE DATE DUE DATE DUE

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

2/05 c:/ClRC/DateDue.indd-p.15

 

PSEUDOELASTIC DEFORMATION OF SHAPE MEMORY
ANNULAR PLATES SUBJECT TO NORMAL TRACTIONS

By

Yuwei Chi

A DISSERTATION

Submitted to
Michigan State University
in partial fulfillment of the requirements
for the degree of

DOCTOR OF PHILOSOPHY

Department of Mechanical Engineering

2005

ABSTRACT
PSEUDOELASTIC DEFORMATION OF SHAPE MEMORY ANNULAR
PLATES SUBJECT TO NORMAL TRACTIONS
By
Yuwei Chi

The pseudoelastic deformation of shape memory annular plates subject to traction
loads is investigated. Pseudoelastic behavior of the shape memory materials iS mod—
eled by a Jg-type deformation theory incorporating an internal variable to accommo-
‘ date the history-dependent hysteresis behavior.

The plane stress and plane strain responses of annular plates subject to forward
transformation are obtained. The effective stresses are found to be non-increasing
along radial direction analytically, which'implies that the maximum effective stress
is always located at the inner boundary of the plate. Therefore the forward transfor-
mation will initiate at the inner boundary and complete at the inner boundary first
upon increasing traction loads. This enables us to study the fully martensite loops
comprised of the load pairs just enough to complete martensitic transformation. Fur-
ther parametric studies on these loops by varying material constants and geometry

are performed. The structure maps partitioning the load in the {1923 p0} plane into

different phase regions are obtained, and the effect of the geometry of the annulus on
the partitioning is studied. The stress fields are Obtained numerically by the shooting
method.

For arbitrary traction histories, an iterative finite difference scheme is developed to
solve the plane stress responses. The stress solution is integrated using the equilibrium
equation and the compatibility equation with given transformation strain field. A
nonlinear equation on transformation strain is Obtained. The equation is discretized
with the finite difference scheme and solved by Newton’s iteration. The numerical
results from the iterations agree with the Shooting results for global loading and global
unloading solutions. The history dependence of stress and martensitic fraction, the
non-monotonousity in effective stress along radial direction and the movement of

phase interfaces are presented.

To my parents

iv

ACKNOWLEDGMENTS

I want to express my sincere gratitude to my major advisor, Dr. Hungyu Tsai, for
his instruction, guidance, and support through the course of my doctoral study and
research. I want to thank my advisor, Dr. Thomas J. Pence, for his support and
constructive suggestions on modeling shape memory materials. I should also extend
my gratitude to Dr. Dahsin Liu, for his valuable advice and help on my research.
I am also grateful to Dr. Zhengfang Zhou, for his valuable discussions on solving
nonlinear ordinary differential equations.

I want to thank Guojing Li, Xu Ding, Yuhui Wang, Chee Kuang Kok, Jun Wu, and
all other colleagues and friends, for their help during my stay at MSU.

Last but not the least, I want to thank my wife Jinghua Xu, who has been taking

great care of me besides working on her own Ph. D. program.

TABLE OF CONTENTS

LIST OF TABLES viii
LIST OF FIGURES ix
1 Introduction 1
1.1 Objectives ................................... 3
1.2 Background .................................. 5
2 Modeling of Pseudoelastic Effects for Shape Memory Materials 16
2.1 Evolution of Martensitic Volume Fraction ................. 17
2.2 Multi—dimensional Stress-Strain Relation .................. 27

3 Plane Stress Analysis of Shape Memory Annular Plates Globally

Loaded from A Phase 31
3.1 Problem Statement .............................. 33
3.2 Reformulation of Governing Equations ................... 37
3.3 A Fully Austenitic Plate (A Phase Distribution) .............. 43
3.4 A Fully Martensitic Plate (M Phase Distribution) ............. 47
3.5 A Plate Containing Mixtures Of Martensite and Austenite (Phase Distri-

butions of type: AC, ACM, CM and C) ................. 52
3.6 Solutions for Global Unloading of M Plates ................ 59
3.7 Discussions .................................. 63

4 Plane Strain Analysis of Shape Memory Annular Plates Globally

Loaded from A Phase 67
4.1 Problem Statement .............................. 68
4.2 Reformulation of Governing Equations ................... 71
4.3 A Fully Austenitic Plate (/1 Phase Distribution) .............. 76
4.4 A Fully Martensitic Plate (M Phase Distribution) ............. 80
4.5 Stress Distributions and Structure Maps .................. 84
5 Plane Stress Analysis of Shape Memory Annular Plates Subject to

Arbitrary Traction History 93
5.1 Problem Statement .............................. 94
5.2 Numerical Solution by Shooting Method .................. 100
5.3 Numerical Solution by Iterative Finite Difference Method ......... 104
5.3.1 Stress Solutions by Integration ...................... 104
5.3.2 Iterative Finite Difference Method .................... 109

vi

5.4 Numerical Results .............................. 120

5.4.1 Comparisons with the Shooting Solution ................. 121
5.4.2 Numerical Results for Plates Subject to Arbitrary Traction History . . 127
5.5 Discussions .................................. 143
6 Conclusions 153
APPENDICES 159
A Mathematica Code for Numerical Examples 159
A1 Structure Map for Plane Stress Analysis .................. 159
A2 Shooting Method for Plane Stress Analysis ................. 160
A3 Structure Map for Plane Strain Analysis .................. 161
AA Shooting Method for Plane Strain Analysis ................. 162
A.5 Iterative Finite Difference Method for Plane Stress Analysis ....... 164
BIBLIOGRAPHY 172

vii

3.1

5.1

5.2

5.3

LIST OF TABLES

The coordinates of the points on Figure 3.13. ...............

Absolute errors and relative errors at t = 1 for the iterating solutions
compared to the shooting solutions ....................

Absolute errors and relative errors at t = 1.5 for the iterating solutions
compared to the shooting solutions ....................

Absolute errors and relative errors for the iterating solutions compared to
the shooting solutions. Time step used is At = 0.05 and totally 101
mesh points (step size h = (r0 — Ti) / 100) are used for simulation.

viii

126

126

145

2.1

2.2

2.3

2.4

2.5

2.6

3.1

3.2

LIST OF FIGURES

(a)Typical one-dimensional stress-strain curve for T > A f. 0571 and of?

are the threshold stresses for the initiation and completion of forward
transformation respectively. a? and O’RT are the threshold stresses
for the initiation and completion of reverse transformation respectively.
(b)Temperature dependence of threshold stresses GET, afT, afifT and

of? When there is no stress, four temperature A f, A3, M3 and

M f mark the transition of forward (A —> M) and reverse (M —> A)
transformation. ..............................

Evolution of martensitic transformation along major loop. ........

Partial forward transformation curves. Regardless of the existing marten-
sitic fraction, partial forward transformation from austenite to marten-
site always starts at 03FT(T) and completes at 0;T(T) .........

Partial reverse transformation curves. Regardless of the existing marten-
sitic fraction, partial reverse transformation from martensite to austen-
ite always starts at 03RT(T) and completes at 0?T(T). ........

Transformation subloops under cyclic loading. ...............

Stabilization of subloops under cyclic loading. ...............

Generic phase distribution within the shape memory alloy annulus. The
inner ring M (r,- S r 3 TM) consists of martensite. The outer ring A
(TA 5 r 3 r0) consists of austenite. The intermediate ring (‘3 (TM <
r < TA) consists of an austenite/martensite mixture. .........

Typical stress distributions (normalized by as) for a plate with M phase
distribution (0(7‘) 2 Of). Here E = 50 GPa, of = 403 = 300 MPa,
ro/ri = 2 and (pl-mo) = (1,5)03. Solid lines represent the case of
a = 0.05. For comparison, stresses for an elastic solution (a = 0)
under the same boundary conditions are plotted as dashed lines. . . .

ix

19

22

25

26

27

27

32

3.3

3.4

3.5

3.6

3.7

3.8

3.9

The FML (Phlly Martensitic Loop) and PAL (Fully Austenitic Loop) for
E = 50 GPa, a = 0.05, of = 403 = 300 MPa and ro/rz- = 2. Tiactions
inside the smaller ellipse (FAL) produce a fully austenitic plate, while
tractions outside the bigger loop (F ML) produce a fully martensitic
plate .....................................

F MLS for various of values. Here E = 50 GPa, a = 0.05, as = 75 MPa
and ro/ri = 2. The smallest loop is the common FAL, and the other
curves are F MLs, each marked with corresponding value of af/os. . .

F MLS for various Ea values. Here a = 403 = 300 MPa and ro/rz- = 2.
The smallest loop is the common AL, and the other curves are FMLS,
each marked with the value of EO/ (500 MPa). .............

FMLS for various ro/rz- values. Here E = 50 GPa, a = 0.05 and 0f =
403 = 300 MPa. The smallest loop is the common FAL, and the other
curves are FMLS marked with ro/ri values ................

Structure map (type I) for ro/rz- = 2. Here E = 50 GPa, a = 0.05 and
of = 403 = 300 MPa. Note that there is no ABM region. Traction
pair (pi, p0) = (5,0.5)03, marked as a dot, produces a BM phase dis-
tribution. The resulting stresses are shown in Figure 3.10. ......

Structure map (type II) for ro/rz- = 5. Here E = 50 GPa, a = 0.05 and
of = 403 = 300 MPa. All six types of phase distributions are present.
Tiaction pair (19,3110) = (5, 0.5)03, marked as a dot, produces a ABM
phase distribution. The resulting stresses are shown in Figure 3.11.

Structure map for the special value ro/ri = 3.58. Lower values produce
type I structure maps, while higher values produces type II. Here E =
50 GPa, a = 0.05 and 0f = 403 = 300 MPa. There is no ABM region,
and ADL contacts MEL at two points. One of the two points (shown
as a dot) represents (pbpo) = (3.565,0.995)03, producing a B phase
distribution with 0(7'2') = of and 0(1‘0) = 03. The resulting stresses
are shown in Figure 3.12 ..........................

3.10 Stresses (normalized by as) for a BM phase distribution. Here E =

50GPa, a = 0.05, of = 403 = 300MPa, ro/ri = 2 and (pl-,po) =
(5, 0.5)03. An elastic solution (a = 0) is shown for comparison (dashed
lines). ...................................

3.11 Stresses (normalized by 03) for a ABM phase distribution. Here E =

50GPa, a = 0.05, of = 403 = 300 MPa, ro/rz- = 5 and (111,130) =
(5, 0.5)03. An elastic solution (a = 0) is shown for comparison (dashed
lines). ...................................

50

51

52

53

55

56

57

58

 

3.12 Stresses (normalized by as) for a Special B phase distribution with TM = r,-

and TA 2 To, for E = 50 GPa, a = 0.05, of = 403 = 300 MPa, ro/ri =
3.58 and (pi,po) = (3.565,0.995)03. An elastic solution (a = 0) is
shown for comparison (dashed lines). ..................

3.13 Correlation of effective stress distributions with a type II structure map.

The parameters used are: E = 50 GPa, a = 0.05, Of = 403 = 300 MPa
and 7‘0/7‘2' = 5. The effective stress distributions corresponding to
various loads in the structure map are shown as subplots. Points 1—10
represent loads on the four bounding loops, marking the transition from
one phase distribution type to another. Typical phase distributions are
shown for loads at points a—f. Values of (pi, 190) used for these points
are listed in Table 3.1 ...........................

3.14 Structure maps for loading and unloading. Here E = 50 GPa, a = 0.05

4.1

4.2

4.3

4.4

4.5

4.6

and of.2T = 00 = 75MPa, a? = 200,05T = 300,0?T = 400, and
r0 = 213. The solid curves correspond to the loading structure map,

and the dashed curves correspond to the unloading structure map. At
point 1, p,- = 190 = 0; at point 2, Pi = 200 and p0 = 400. .......

The FML (Fully Martensitic Loop) and PAL (Fully Austenitic Loop) for
E = 50 GPa, a = 0.05, V = 0.3, of = 403 = 300 MPa and ro/rz- = 2.
Tractions inside the smaller ellipse (FAL) produce a fully austenitic
plate, while tractions outside the bigger loop (FML) produce a fully
martensitic plate. Part (b) is the enhanced view of the zone enclosed
in the dashed rectangle in (a). ......................

FMLS for various of values. Here E = 50 GPa, a = 0.05, 1/ = 0.3, as =
75 MPa and r0 /1‘,- = 2. All the FMLS are marked with corresponding
value of af/os. ..............................

FMLS for various Ea values. Here V = 0.3, of = 403 = 300 MPa and

ro/ri = 2. All the FMLS are marked with the value of EO/ (500 MPa).

FMLS for various V values. Here E = 50 GPa, a = 0.05, Of = 403 =

300 MPa and ro/rz- = 2. All the FMLS are marked with the value of V.

FMLS for various ro/rz- values. Here E = 50 GPa, a = 0.05, V = 0.3 and
of = 403 = 300 MPa. All the FMLS are marked with ro/ri values. . .

Structure map (type I) for ro/ri = 2. Here E = 50 GPa, a = 0.05, l/ = 0.3
and of = 403 = 300 MPa. Part (b) is the enhanced view of the zone
enclosed in the dashed rectangle.in (a). Note that there is no ABM
region. Traction pair (pl-mo) = (8,1)03, marked as a dot, produces a

BM phase distribution. The resulting stresses are shown in Figure 4.9.

60

65

66

82

83

84

85

86

87

4.7 Structure map (type II) for ro/ri = 8. Here E = 50 GPa, a = 0.05 ,
V = 0.3 and 0f = 403 = 300 MPa. Part (b) is the enhanced view of the
zone enclosed in the dashed rectangle in (a); part (c) is the enhanced
view of the zone enclosed in the dashed rectangle in (b). All six types
of phase distributions are present. Traction pair (piano) = (8,1)03,
marked as a dot, produces a ABM phase distribution. The resulting
stresses are shown in Figure 4.10. ....................

4.8 Structure map for the special value I‘D/r, = 5.52. Here E = 50 GPa,
a = 0.05, V = 0.3 and of = 403 = 300 MPa. Lower ro/rz- values
produce type I structure maps, while higher ro/rz- values produces type
II. Part (b) is the enhanced view of the zone enclosed in the dashed
rectangle in (a); part (c) is the enhanced view of the zone enclosed
in the dashed rectangle in (b). There is no ABM region, and ADL
contacts MEL at two points. One of the two points (shown as a dot)
represents (pi,po) = (4.43,0.68)03, producing a B phase distribution
with C(73) = 0f and 0(r0) = as. The resulting stresses are shown in
Figure 4.11. ................................

4.9 Stresses (normalized by as) for a BM phase distribution. Here E =
50 GPa, a = 0.05, V = 0.3, of = 403 = 300MPa, ro/ri = 2 and
(pi, p0) = (8,1)03. An elastic solution (a = 0) is Shown for compari-
son (dashed lines). ............................

4.10 Stresses (normalized by as) for a ABM phase distribution. Here E =
50GPa, a = 0.05, V = 0.3, (If = 403 = 300MPa, ro/rz- = 8 and
(11,3190) = (8, 1)os. An elastic solution (a = 0) is shown for comparison

(dashed lines) ................................

4.11 Stresses (normalized by as) for a special B phase distribution with TM = r,
and 7'A = To, for E = 50 GPa, a = 0.05, V = 0.3, of = 405 = 300 MPa,
ro/ri = 5.52 and (15,120) = (4.43,0.68)03. An elastic solution (a = 0)
is shown for comparison (dashed lines) ..................

5.1 Load steps in (19,3120) plane. Here E = 50 GPa, a = 0.05, aftT = 00 =
100MPa, a? = 200, 057‘ = 300, of? = 400, and r0 = 213. The
curves on the plate are from the loading structure map. A total of 40
time steps (At = 0.05) is used. For 0 < t S 1 (step 1 to step 20), the

plate is globally loaded to ABM phase combination. For 1 < t S 2
(step 21 to step 40), the annulus is globally unloaded to A phase.

xii

88

90

91

91

92

101

5.2

5.3

5.4

5.5

5.6

5.7

Effective stresses (normalized by 00) by the shooting method for the load-
ing steps corresponding to Figure 5.1. A total of 40 time ste is used,
and E: 50 GPa, a — _0.,05 of.” = 00— — 100MPa, as =200,
JET = 300, of? = 400, and r0 = 21",. (a) Effective stresses for

0 S t S 1 (initial step to step 20), the plate undergoes global loading.

(b) Effective stresses for 1 _<_ t S 2 (step 20 to step 40), the plate

undergoes global loading. ........................

Evolution of martensitic fractions for r = 1, r = 1.5, and 1' = 2. Effective
stress a is normalized by 00. The applied load steps correspond to

F Igure 5.1. A total of 40 time step is usedde E = 50 GPa, a --= 0.05,
of — 00:100MPa,o£:"T = 200, of = 300, of?“ = 400, and
r0— —— 2r,- ...................................

Comparison Of stress distributions (normalized by 00) between the iter-
ating solution and the shooting solution at t= 1. Here E = 50 GPa,
a— _ 0. 05, ofT = ao_=100MPa,a§T= 200,0 aFT = 300,011.FT = 400,
and r0— — 273. The applied tractions correspond to t = 1 are p,- = 500
and p0 = 300. The solid curves correspond to the shooting solution,
while the scatter points are results from the iterating solution with
mesh points of 5, 9, and 17 respectively. ................

Comparison of stress distributions (normalized by 00) between the iter-
ating solution and theT the shooting solution at t- -— 1.5 Here E =
50GPa, a— — 0.05, afT = 00— = 100MPa, agiT— = 200, ofT = 300,

afT— — 400, and r0— — 27“,. The applied tractions correspond to t- — 1. 5
are p,- = 2.500 and p0 = 1.500. The solid curves correspond to the
shooting solution, while the scatter points are results from the iterating

solution with mesh points of 5, 9, and 17 points respectively. .....

Traction history of Of Example 1. The load follows a path 0 —+ a —-» b ->
c —> a in (19,3190) plane. Other thin curves are part of the loading
structure map. This load history will be referred to as path I.

Comparison of effective stress distributions (normalized by 00) between
t = 1 and t = 4 for to traction path I. The applied tractions at t = 1
and t = 4 are both p,- = p0 = 200. The scatter points and solid curve
are effective stresses at t = 4. The dashed-line is the effective stress at
t = 1 and 0(r) = 200 ............................

103

123

124

129

130

5.8

5.9

5.10

5.11

5.12

5.13

5.14

5.15

5.16

Comparison of radial stress and hoop stress distributions (normalized by
00) between t = 1 and t = 4 for traction path I. The applied tractions
at t = 1 and t = 4 are the same, with p,- = p0 = 200. The scatter
points and solid curves are stresses at t = 4. The dashed-line is the
stress at t = 1 and (rm-(r) = 099(7‘) = —200. ..............

Time history plot of effective stresses (normalized by 00) at r = 1, 7' = 1.5
and r = 2. Applied traction is path I, and a time step At = 0.05 and
101 mesh points are used in the computation. .............

Evolution of martensitic fractions at r = 1, r = 1.5 and r = 2. Effective
stress a is normalized by 00. Applied traction is path I, and a time
step At = 0.05 and 101 mesh points are used in computation.

Effective stresses (normalized by 00) for 2 S t S 3. Applied traction
is path I, and a time step At = 0.05 and 101 mesh points are used
in computation. During this period, no forward transformation takes
place. Part of the plate undergoes reverse transformation for 2.85 S
t < 3 .....................................

Difference in effective stresses (normalized by 00) between two consecutive
steps (Ao[t] = a[t] — o[t — At]) for the period 2 < t S 3. Applied
traction is path I, and a time step At = 0.05 and 101 mesh points are
used in computation. ..........................

Difference in martensitic fractions between consecutive steps (A5 [t] =
6 [t] - 5 [t — At]) for the period 2 < t S 3. Applied traction is path
I, and a time step At = 0.05 and 101 mesh points are used in compu-
tation. ...................................

Effective stresses (normalized by 00) for the period 3 S t S 4. Applied
traction is path I, and a time step At = 0.05 and 101 mesh points are
used in computation. ..........................

Difference in effective stresses (normalized by 00) between consecutive
steps (A0[t] = o[t] — a[t - At]) for the period 3 < t S 4. Applied
traction is path I, and a time step At = 0.05 and 101 mesh points are
used in computation. ..........................

Difference in martensitic fractions between consecutive steps (A5 [t] =
6 [t] - g [t — At]) for the period 3 < t S 4. Applied traction is path
I, and a time step At = 0.05 and 101 mesh points are used in compu-
tation. ..................................

133

134

5.17 Traction history for Example 2. This traction history is denoted as path H.140

xiv

5.18 Time history plot of effective stresses (normalized by 00) at r = 1, r = 1.5,
and r = 2. Applied traction is path II, and a time step At = 0.01 and
101 mesh points are used in computation. ............... 140

5.19 Traction history for Example 3. This traction history is denoted as path III. 141

5.20 Time history plot of effective stresses (normalized by 00) at 1‘ = 1.0, 1.2, . . . ,2.0.
Applied traction is path III, and a time step At = 0.01 and 201 mesh

points are used in computation. ..................... 146
5.21 Evolution of martensitic fractions at r = 1.0, 1.2, . . . ,2.0. Effective stress

a is normalized by 00. Applied traction is path III, and a time step

At = 0.01 and 201 mesh points are used in computation. ....... 147

5.22 Convergence of effective stresses (normalized by 00). Applied traction is
path III, and a time step At = 0.01 and 201 mesh points are used in
computation. (a) Effective stresses at t = 1,3, . . . , 19. (b) Effective
stresses at t = 2, ,...,20. ........................ 148

5.23 Movement of r A (B /A interface) and TM (M/ B interface). Applied trac-
tion is path V, and a time step At = 0.01 and 201 mesh points are
used in computation. ........................... 148

5.24 Traction history for Example 4. This traction history is denoted as path V. 149

5.25 Time history plot of effective stresses (normalized by 00) at r‘ = 1.0, 1.2, . . . , 2.0.
Applied traction is path V, and a time step At = 0.01 and 201 mesh

points are used in computation. ..................... 150
5.26 Evolution of martensitic fractions at r = 1.0, 1.2, . . . ,2.0. Effective stress

a is normalized by 00. Applied traction is path V, and a time step

At = 0.01 and 201 mesh points are used in computation. ....... 151

5.27 Movement of TA (B /A interface) and TM (M/ B interface). Applied trac-
tion is path V, and a time step At = 0.01 and 201 mesh points are
used in computation. ........................... 152

XV

Chapter 1

Introduction

Shape memory materials refer to the class of materials that exhibit shape memory
efiect and pseudoelasticity. Within certain temperature range, if a shape memory ma-
terial initially undergoes deformation due to mechanical load, and the deformation
does not recover completely after the applied load is removed, it can return to the
previous undeformed shape by heating. This is called shape memory effect. At tem-
peratures higher than a critical value (depending on the material), a shape memory
material can have a large strain (usually up to 5—10%) upon loading, and the strain
recovers completely upon unloading. However, the loading stress-strain curve and the
unloading stress-strain curve follow different paths. This property of shape memory
materials is called pseudoelasticity (or superelasticity). The difference in loading and
unloading curves distinguishes pseudoelastic behavior from general nonlinear elastic
behavior, in which the loading and unloading stress-strain relations are identical.

The underlying mechanism of shape memory effect and pseudoelasticity is the diffu-

sionless martensitic phase transformation, in which austenite phase (A phase) con-

verts to martensite phase (M phase) or vice versa. The .A phase favors high temper—
ature and low stress, conversely the M phase favors low temperature and high stress.
A detailed description of the macroscopic behavior of shape memory materials and
the crystallographical representation of martensitic transformation can be found in
Otsuka and Wayman (1998).

Thermoelastic martensitic transformation may take place when shape memory mate-
rials are heated up or cooled down. The transformation is reversible and it is denoted
as temperature-induced martensitic transformation. If no stress is applied to shape
memory materials, the temperature-induced martensitic transformation is described
with the aid of the four special transformation temperatures Ms, M f, A f and A3.
The austenite phase is converted to the martensite phase by cooling the material
through the temperature interval MS > T > M f (the subscripts refer to “start” and
“finish”), and it is denoted as forward transformation. The material is converted
back to the austenite phase by heating the material through the temperature interval
A3 < T < A f, and it is denoted as reverse transformation. Typically M3 < A3 so
that M f < M5 < A,» < A f. The difference in these transformation threshold tem-
peratures creates hysteresis between forward martensitic transformation and reverse
martensitic transformation due to temperature change.

Similarly, stress-induced martensitic transformation may take place if external loads
are applied to shape memory materials. In uniaxial extension experiment of shape
memory wires under isothermal condition with a temperature higher than A f, four
threshold stresses, 0.5T, affrT, 0.5T and of]: mark the starting and finishing of

forward transformation and reverse transformation. The superscript FT refers to

2

forward transformation, and the superscript RT refers to reverse transformation.
When the stress applied to the materials increases, the forward transformation starts

at GET and completes at ofT; when the stress applied to the materials decreases, the

reverse transformation starts at a? and completes at 0?? Typically 0§T < JET

FT

RT<0§T<US

hence a f < of? Hysteresis also exists between the forward
and reverse transformations due to change in stress. Pseudoelasticity is the result of

stress-induced martensitic transformation.

1 .1 Objectives

The purpose of this work is to study some boundary value problems of shape memory
materials. Specifically, we are interested in how annular plates respond to normal
tractions applied on their edges.

This research follows closely recent structural analyses of ring structures associ-
ated with martensitic transformation. Giannakopoulos and Olson (1993) investi-
gated the axisymmetric deformation of ring structures for materials that admit dila-
tional martensitic transformation instead of the shearing martensitic transformation
in shape memory materials. Dilational martensitic transformation is responsible for
transformation toughening in certain ceramic materials (Budiansky, et al., 1983). The
martensitic transformations associated with shape memory behavior in metal alloys
have little, if any, volume change.

Birman (1999) considered an infinite plate with a traction free hole under stress-

induced transformation. The transformation strain is assumed to be isotropic expan—

sion, and only forward transformation under plane stress setting is considered in the
analysis. Part of the analysis uses an additional assumption that the ratio of radial
stress to hoop stress is the same as that of the elastic solution, which is questionable
as the transformation strain field may change the distribution of the stresses. Also the
geometry and boundary conditions used in the analysis are special, and the results
may not be applicable to finite size ring structures.

Here we analyze the pseudoelastic deformations for shape memory annular plates
using a relatively more refined material model, as compared to the one used by Bir-
man (1999). A pseudoelastic model capable of describing the transformation behavior
under various stress load history is used, hence the martensitic transformation is not
limited to forward transformation. The transformation strain is described by J2 de-
formation theory hence it is volume preserving, and this is consistent with the shear
nature exhibited in martensitic transformation. The traction loads are applied on
both the inner and outer edges of the plates and can have any combination. First,
the plane stress and plane strain analyses for plates undergoing forward transforma-
tion from austenite phase are studied, and parametric studies on the effect of material
constants and plate geometries on the responses of the plates are performed. Then
the plane stress analyses of annular plates under various traction histories, whereby
both forward and reverse transformation may take places within the plate, are carried

out.

 

1.2 Background

Since the 19508, shape memory materials have drawn close attention primarily for
their memory effect and pseudoelasticity, which potentially have a wide range of
applications (Wayman, 1980). For the past decades, shape memory materials are
increasingly used as active materials, energy dissipation materials, and biomedical
materials (Otsuka and Wayman, 1998). They have been very successfully applied
to active or passive vibration control of civil structures (Saadat et al., 2002), air-
conditioner flaps, microactuators and sensors (Otsuka et al., 2002), and medical im-
plants such as stents (Duerig, 2002). The design of shape memory devices is based
on either shape memory effect or pseudoelasticity, with devices relying on pseudoe-
lasticity being more successful than those on shape memory effect (Otsuka, 2002).
It is essential to understand how the devices respond to temperature changes and
mechanical loads in design. This research on shape memory annular plates may be
useful in shape memory applications involving ring or tube structures.

The early studies on shape memory materials are more of experimental investiga-
tions. Most of the experiments focused on identifying materials exhibiting shape
memory effect and pseudoelasticity, observing microstructure change during marten-
sitic transformation, or calibrating macroscopic one-dimensional behavior under var-
ious thermal-mechanical loadings (Delaey et al., 1974).

The threshold temperatures associated with martensitic transformation are important
parameters for shape memory materials. Cornelis and Wayman (1976) measured the

resistance of CuZn alloy subject to various temperature cycles. They found that the

 

threshold temperatures under which the forward transformation or the reverse trans-
formation completes are constant regardless of the degree of previous transformation.
The evolution of martensitic volume fraction due to temperature change is enclosed
in an envelope formed by two paths; one path corresponds to forward transformation
starting from austenite, the other path corresponds to reverse transformation starting
from martensite. We will use transformation major-loop to denote these two paths,
and transformation subloop to denote other paths inside the envelope in the sequel.
In addition to the invariance in transformation-completing threshold temperatures
for subloops, Paskal and Monasevich (1981) demonstrated that the threshold tem-
peratures for initiating the forward transformation or the reverse transformation are
also independent of the degree of previous transformation in TiNi alloy. Those results
are widely adopted in modeling of temperature-induced transformation behavior.

As the counterpart of the threshold temperatures in temperature-induced transfor-
mation, the threshold stresses associated with pseudoelastic deformation were also
studied. Tobushi et al. (1993) performed tensile tests on TiNi shape memory wires
with various maximum strains. The threshold stresses for forward and reverse trans-
formation are found to be linearly related to temperature. The elastic modulus for
the phase mixture can be approximated by an averaging formula using phase frac-
tions with the assumption that the same stress is applied on both phases. Tanaka
et al. (1995a) obtained similar results from experiments on TiNi polycrystalline and
Fe-based polycrystalline alloys. The transformation threshold stresses for starting
and finishing of martensitic transformation are linearly dependent on temperature

for both alloys. For TiNi alloy, the stress ranges for both forward transformation and

6

reverse transformation are fairly small. They confirmed that the threshold tempera-
tures for TiNi alloys are essentially independent of previous thermal loading history,
while pointing out that the threshold temperatures for Fe-based alloys are strongly
history dependent.

Novak et al. (1999) performed tension/ compression tests on Cu-Al-Ni single crys-
tal shape memory under different temperatures. The alloys exhibit anisotropy and
asymmetry in tension / compression, and the shape of hysteresis loops is strongly tem-
perature dependent. The tension-torsion tests on thin-wall tubes by Lim and Mc-
Dowell (1999) revealed the tension/compresion asymmetry and path-dependent be-
havior in martensitic transformation. The nonproportional responses of the material
are found to be very complicated, as the martensitic transformation and martensitic
reorientation may take place at the same time. Developing multiple-dimensional con-
stitutive models to represent the complicated stress-strain relation in shape memory
materials is challenging.

Quite a few models were developed to account for the behavior of shape memory
materials exhibited in experiments. Martensitic transformation associated with shape
memory effect and pseudoelasticity are often modeled with the aid of free energy.
Falk (1980) proposed a free energy as a polynomial function of shear strain and
temperature based on the Landau theory. Stress-induced and temperature—induced
transformations and shape memory effect can be derived from the minimization of
this specific energy. The results match experimental data qualitatively.

Muller and Xu (1991) used interfacial energy to interpret the pseudoelastic behavior.

The total free energy is taken as the sum of the interfacial energy and individual

7

free energy of the both phases weighted by individual phase volume fraction. The
individual free energy is specialized from the one by Falk (1980) and defined as a
quadratic function of strain and temperature, and the interfacial energy is linear to
the product of the two phase fractions. Minimization of the total free energy gives
the stress-strain hysteresis loops. The width of the hysteresis is determined by the
interfacial energy, and it grows slightly as temperature increases. The interior loops
associated with cyclic loading is accommodated by introducing metastable states,
and phase transformation occurs when sub100ps cross the unstable phase equilibrium
line. Although the results seem to agree with the stress-strain curves from their
deformation-controlled experiment on CuZnAl alloy, Tanaka et al. (1995a) questioned
the existence of such a single unstable line based on their experimental results.
Kamita and Matsuzaki (1998) studied one-dimensional pseudoelastic theory using a
similar free energy by Muller and Xu (1991), with interfacial energy changed to in-
teraction energy. The interaction energy is determined by martensitic fraction and
is different between forward transformation and reverse transformation. The total
free energy is used to construct extended potential energy for deriving the equilib-
rium condition of the two-phase state. The results appear to be consistent with the
experimental results of Tobushi et al. (1993).

Bernardini and Pence (2002a) developed a thermodynamic model involving both a
free energy function and a dissipation function. The free energy function is different
from those used in the above models. A common free energy is used for both austenite
and martensite phases, but different dissipation functions are adopted. The constitu-

tive equations are determined from these free energy and dissipation functions based

8

on laws of thermodynamics. The model is able to capture temperature-induced phase
transformation and stress-induced phase transformation under different thermal con-
ditions.

Without resorting to free energy functions, phenomenological description of macro-
scopic transformation behavior is straightforward and easy to use. Tanaka (1990)
described the thermomechanical behavior of shape memory materials phenomeno-
logically by giving martensitic fraction evolution functions for forward and reverse
transformation explicitly. Liang and Rogers (1990) improved the model to include
the effects of initial phase mixture, and replaced the exponential functions with the
cosine functions to describe the martensitic fraction evolution during transformation.
Ivshin and Pence (1994a) proposed a rate-independent model to study the transforma-
tion subloops and hysteresis. For forward transformation, the martensitic volume frac-
tion accumulates at a rate proportional to the austenitic volume fraction; for reverse
transformation, the austenitic volume fraction accumulates at a rate proportional
to the martensitic volume fraction. Transformation starting and finishing threshold
temperatures are taken as constants. A general phase evolution function based on the
Duhem—Madelung framework is proposed to describe the rate-independent transfor-
mation behavior. The conditions for the existence and uniqueness of transformation
subloops when temperature alternates between two fixed values are obtained. In
addition, restricting conditions on the model to preserve the monotonicity in trans-
formation and containment of subloops within major loops are presented in detail.
These issues are important for a well-defined hysteresis model. The drifting of phase

fraction as a result of oscillatory temperature cycles and geometrical features associ-

9

 

ated with subloops are also discussed. The framework provided by this work is quite
general and can be easily used to derive specific rate-independent models. This model
is used in our analysis with minor modifications.

The above model was extended by Ivshin and Pence (1994b) to accomodate phase
evolutions under temperature changes and stress loads. The temperature and stress
are combined into one driving input for the phase transformation. The phase evolu-
tion functions preserve similar formulae as those in Ivshin and Pence (1994a), with
temperature changed to this new driving input variable. The model by Liang and
Rogers (1990) is a special case of this model. Stress-strain relations associated with
different thermal conditions and mechanical loads are presented. This single-variant
model is further extended to a two-variant model by Wu and Pence (1998). By unfold-
ing a phase diagram triple point on the stress-temperature plane, general phase dia-
grams describing the conversion between any two of the three compositions (austenite
and two martensite variants) are presented. Two driving input functions, both func-
tions of temperature and stress, are introduced to describe the evolution of phase
fractions under various temperature and stress histories. This model is able to ac-
count for the detwinning behavior, which single-variant models are unable to capture.
Several multi-dimensional models have been proposed. Many models take a form
modified from plasticity. Due to the complexity in transformation behavior, the con-
stitutive relations are often expressed in rate form. However, analytical investigation
on an incremental constitutive model is, in general, very difficult.

Bondaryev and Wayman (1988) presented a three-dimensional model to describe pseu-

doelastic deformation and shape memory effect. At high temperature, the Gibbs free

10

energy of either phase equals the summation of strain energy and chemical energy,
and the latter is a function of temperature. The change in the free energy during
transformation is simplified by neglecting hydrostatic terms in the strain energy and
linearizing the chemical energy about temperature. The condition for initiating phase
transformation is obtained with the assumption that such an energy change equals a
constant reflecting the resistance in transformation. Temperature dependent surfaces
in stress space are used to determine the transformation strain accumulation and re-
orientation, and the transformation strain increments follow J2 flow theory. At low
temperature, ferroelasticity is attributed to the reorientation of martensite plates.
Numerical results on stress-strain-temperature relations are given for the biaxial load
case. Guided by this work, we use a J2 deformation theory to describe the transfor-
mation strain, and extend the one-dimensional transformation behavior described by
Ivshin and Pence (1994a) to a three-dimensional setting.

Gillet et al. (1995) proposed a three-dimensional model that reflects the tension-
compression asymmetry by introducing a yield function dependent on J2, J3, trans-
formation strain and temperature. The introduction of J3 is important for the purpose
of addressing asymmetry in tension-compression, as compared to the conventional J2
theory. The transformation strain increment is obtained using the associated rule, and
hysteresis is introduced by providing different loading and unloading yield functions.
Boyd and Lagoudas (1996) developed a three-dimensional model using Gibbs free
energy and dissipation potential. The total free energy equals the summation of indi-
vidual free energy of different species and mixing free energy in the crystal. The dis-
sipation potential is approximated by a quadratic function of the net thermodynamic

11

force associated with entropy production, and is used for deriving the transformation
kinetic equations and the reorientation strain rate. The total transformation strain
rate combines the reorientation part and transformation part, with the latter coaxial
to a modified stress deviator tensor. The model involves quite a few constants to be
determined experimentally, however, the author didn’t provide a comparison of the
model with any multi-dimensional experimental result. Qidwai and Lagoudas (2000)
derived a J2 — J3 — 11 model based on the concept of maximum transformation dissi-
pation rate, and the tension / compression asymmetry can be captured by this model.
Numerical examples using this model are given, but are limited to the homogeneous
stress situation. The non-proportional loading cannot be captured by this model as
additional internal state variables may be needed to account for the reorientation of
transformation strain. This is also pointed out in Lim and McDowell (1999). They
demonstrated the discrepancy between the numerical results from a J2 — J3 type
model and the experimental results, and suggested that the modified plasticity model
is insufficient in describing non-proportional loading behavior.

Micromechanics has been used to develop constitutive models for shape memory ma-
terials. Sun and Hwang (1993a, 1993b) proposed a micromechanical model for poly-
crystalline shape memory alloys under thermal or mechanical loading. By assuming
transformed grain to be a spheroid, the elastic energy stored in the polycrystalline
is obtained by the Mori-Tanaka self—consistent theory. The effect of microstructural
mechanisms including internal stress, twin boundary and interface friction are ac-
counted forby including surface energy and transformation dissipation energy. Pseu-

doelasticity and shape memory effect due to transformation and reorientation are

12

studied, and temperature, stress and loading history-dependent yielding surfaces are
given. Reorientation and pseudoelastic cycles under non-proportional loading paths
are also discussed.

Huang and Brinson (1998) developed a multi-variant model for single-crystal shape
memory materials. The material is treated as an austenite matrix with randomly
distributed martensite inclusions. These inclusions are groupings of martensitic vari-
ants free of interaction energy within the group. The interaction energy is computed
by the Eshelby-Kroner approach. This model is able to capture temperature-induced
transformation and certain tension/ compression asymmetry.

Lim and McDowell (2002) used the finite element method to implement a microme-
chanical model for textured shape memory polycrystals. Intergranular interactions
are simulated by choosing appropriate finite element mesh, and user material subrou-
tine is developed to specify the constitutive relations. The model is able to capture
tension/ compression asymmetry under the uniaxial setting and the rate-dependence
due to latent heat during transformation. Although that these micromechanical mod-
els seem to represent the experimental results better than other three-dimensional
models described above, it is impractical to use these models in our analysis where
an explicit stress-strain relation is desired.

The above review of publications in the modeling of shape memory materials is by
no means complete, as only those closely related to our work are included. More

detailed discussions on modeling considerations and comparison of different models

can be found in Otsuka (1986), Birman (1997), and Bernardini and Pence (2002b).

13

Far from the richness of publications on modeling, only a few structural analyses on
shape memory materials were reported. Brinson and Lammering (1993) deveIOped a
one-dimensional finite element algorithm for shape memory alloys. A pipe structure
with an inner ring composed of normal elastic material, and an outer ring composed
of shape memory material, is analyzed where the outer ring is approximated by one-
dimensional truss elements. Such a simplification is quite crude, and if the ratio of the
outer radius to the inner radius of the outer ring is large, a cluster of truss elements
is not a good representation of this ring structure.

Kuczma and Mielke (2000) used the one-dimensional model by Huo and Muller (1993)
to analyze a bar with non-uniform cross-sectional area subject to displacement bound-
ary conditions. Note for this analysis, although the stress is also non-uniform along
the bar, it is inversely proportional to the cross-sectional area. Therefore, knowledge
of the stress at any location in the bar is sufficient to determine the stress everywhere
in the bar.

Berman and White (1996) studied the responses of a shape memory composite under
thermoelastic transformation. The strains due to temperature-induced transforma-
tion and thermal expansion are combined to a thermal strain. The relation between
this thermal strain and temperature are obtained from the experiment on shape mem-
ory wires. Therefore, it is essentially using an iso—expansion model for transformation
strain while neglecting the thermal expansion, which is much smaller than the trans-
formation strain. Neither stress-induced transformation nor partial transformation

following subloops are considered in the analysis.

14

Trochu and Qian (1997) used bilinear loading and unloading functions to describe
martensitic transformation, and transformation strain is obtained using plastic flow
theory. Two numerical examples, bending of a beam and compression of a spring
disk, are given. Assuming that the transformation strain is isotropic expansion and
proportional to the first invariant of stress tensor, Yan et al. (2002) performed the
finite element analysis on pseudoelastic deformation for a crack with ABAQUS. Only
major transformation loops are considered. The numerical results were used to per-
form toughness analysis. Although complicated models can be implemented with
finite element method, it only gives the responses of the structure under specific load-
ing conditions. In general, it is hard to perform parametric studies on the effects of

the material properties and geometry based on results from finite element method.

15

Chapter 2

Modeling of Pseudoelastic Effects

for Shape Memory Materials

In this chapter, we describe the constitutive model for pseudoelastic deformation
in shape memory materials adopted for future analysis. The one-dimensional stress-
strain relation and the evolution of martensitic transformation for the one-dimensional
setting are given, and they are extended to the three-dimensional setting based on
the J2-deformation theory.

Martensite has a microstructure consisting of a number of crystallographic variants
that are related to each other by various crystallographic rotations and reflections. At
the engineering scale, the transformation strain can be approximated by a function
of particular variant collections, and here we use the martensitic volume fraction,
denoted by E, as an internal variable to describe the collection of martensite variants.
For the total strain during phase transformation, a standard treatment is to decom-

pose the strain tensor 5 into elastic strain eel and transformation strain E”, with

16

elastic strain eel observing Hooke’s law and transformation strain a” determined by
martensitic phase fraction and stresses.

For the one-dimensional setting, the total strain can be written as

l

e = 58 + a" with eel = % and a" = a5, (2.1)

where 5 is the axial total strain, 561 the axial elastic strain, a" the axial transformation
strain, a the axial tensile stress, E the Young’s modulus, and a the material constant
corresponding to the maximum transformation strain in uniaxial tension tests. For
this simplified model, the Young’s moduli of martensite and austenite are assumed to
be equal. A refinement of this model can be made by averaging the moduli of both
phases according to their volume fractions (Tobushi et al., 1993).

How the martensitic volume fraction 6 evolves in response to stress and temperature
change is essential in this one-dimensional model. The evolution of the martensitic

fraction under uniaxial load is given in the next section.

2.1 Evolution of Martensitic Volume Fraction

The martensitic phase fraction E is usually a function of stress and temperature, and
is history dependent. As mentioned in the introduction, four special transformation
temperatures Ms, M f, A f and A3 are associated with the starting and finishing
temperatures of stress-free martensitic transformation. At temperatures above the

austenite finish temperature A f, similarly, four threshold stresses JET, ofT, ogT

17

and 0?T are associated with stress-induced phase transformation. Pseudoelastic

response in general refers to the isothermal behavior of shape memory materials at
temperatures above A f due to the generation of stress-induced martensite. This
process has been extensively studied and modeled in one-dimensional settings such
as uniaxial tension of thin wires.

The stress-strain curve for increasing stress in uniaxial extension is characterized by
a pseudoelastic loading plateau as the material converts from austenite to martensite
(the forward transformation). A schematic plot of stress-strain curve under a uniaxial
extension test is shown as the upper curve in Figure 2.1(a). The horizontal grid
lines mark the starting and finishing stresses of forward transformation and reverse
transformation. This conversion takes place over a relatively small stress interval
from 05T(T) to ofT(T) which depends on temperature (whence the argument T).
The superscript FT denotes forward transformation. The threshold stresses GET and
of], increase with temperature at an essentially constant rate so as to be consistent
with the Clausius—Clapeyron relation. This linearity was observed in experimental
investigations, such as those by Tobushi et al. (1993) and Tanaka et al. (1995a).
Figure 2.1(b) shows a schematic plot of threshold stresses as functions of temperature.
For the description of phase transition at temperature below A f where shape memory
effect exists, please refer to Bernardini and Pence (2002b).

Similarly, the stress-strain curve for decreasing stress in uniaxial tension at T > A f
is characterized by a pseudoelastic unloading plateau (below the loading plateau)
as the material converts martensite back to austenite (the reverse transformation),

which is shown as the lower curve in Figure 2.1(a). Threshold stresses are 033T and

18

 

0'

 

 

 

 

6}” at?" oi? 6?
(a) (13)
FT FT

Figure 2.1: (a)Typical one-dimensional stress-strain curve for T > A f. 03 and a f

are the threshold stresses for the initiation and completion of forward transformation
respectively. of” and 011,271 are the threshold stresses for the initiation and comple-

tion of reverse transformation respectively. (b)Temperature dependence of threshold

stresses GET, OFT, of” and URT. When there is no stress, four temperature A f,

A3, M3 and M f mark the transition of forward (A —> M) and reverse (M —) A)
transformation.

0?, where again the temperature dependence is essentially linear. The superscript
RT denotes reverse transformation. In particular, the equality, of.” (A f) = 0, re-
flects the fact that unloading at temperatures below A f results in an incomplete
reverse transformation. The remaining threshold stress functions 033T (T), 03FT(T)
and of? (T) are defined on the respective temperature intervals T 2 A3, T 2 M3,
T 2 Mf. These functions obey gyms) = 0, 05T(Ms) = 0 and ofT(Mf) = 0.
The temperature ordering M f < M3 < A; < A f implies that for T > A f we have

ame < asRTm < achr) < achr).

19

In Chapter 1, we briefly defined transformation major loop and transformation subloop
associated with temperature—induced transformation. These two terms are now ex-
tended to stress-induced transformation. The evolution of martensitic fraction about
stress on a {-0 plot is enclosed in an envelope, which is the transformation major
loop. The enclosed paths are transformation subloops. For conciseness, we omit the
word “transformation” and simply use major loop and subloop to describe these two
types of paths. The martensitic transformation along the major loop is called full
transformation, and the martensitic transformation along subloops is called partial
transformation. The containment of subloops inside the major loop is consistent
with one-dimensional experiments, and it is regarded as an important condition in
modeling martensitic transformation (Ivshin and Pence, 1994).

The envelope can be defined by two functions gF T(o, T) and gRT(o, T) Different
schemes for subloops can be developed to refine the hysteresis behavior of martensitic
transformation. There are several models that considered the hysteresis behavior
involving subloops, for example, Liang and Rogers (1990), Brinson and Lammer-
ing (1993), and Ivshin and Pence (1994a). The treatment by Ivshin and Pence (1994a)
is more rigorous than the other two, and we will follow the framework of this work
to describe the hysteresis behavior.

The envelope function gF T(o, T) corresponds to the forward transformation that
starts from a state of pure austenite (5 = 0) and zero stress. The martensite phase

fraction 5 for forward transformation is

5 = 9“"(0. T). (2.2)

20

Similarly, the envelope function gRT(0, T) corresponds to the reverse transformation
that starts from a state of pure martensite (E = 1). The martensite phase fraction 6

for reverse transformation is

5 = 9RT(0,T)- (2-3)

A typical example of gFT(o, T) and gRT(o, T) takes the following form (Liang and

Rogers, 1990):

 

 

0 for o S U§T(T),
9FT(U,T) = _1—{1 cos [ o — (jg-TC” al} for 05T(T) < o < 0FT(T),

2 afTa‘) — a§T(T) f

1 for o 2 0?T(T)

(2.4)

and

1 for a Z 0§T(T),

RT
RT T = ‘7 “ ”f (T) RT RT

9 (0, ) §{1— cos [0§T(T) _ 0?T(T)7r]} for of (T) < o < 03 (T),

0 for o S 0?T(T).

(2.5)
Figure 2.2 is a schematic plot for the above functions. These envelope functions only
give the evolution of martensitic phase fraction for forward transformation starting
from austenite phase (5 = 0), or reverse transformation starting from martensite
phase (5 = 1). This is sufficient for the description of complete transformation be-

havior, but inadequate for the description of partial transformation behavior asso-

21

0.5 ‘ gRT(6,T)

 

 

0'

 

T > T
6}”(T) 6? (T) of "(1‘) 6}” (T)
Figure 2.2: Evolution of martensitic transformation along major loop.

ciated with subloops. Extensions that allow for the description of subloops due to
partial transformation followed by transformation residual can be accomplished by
extending the above framework using history dependent kinetic rules that determine
6 obeying gF T(o, T) S .5 S gRT(a, T). One such extension can be found in Ivshin
and Pence (1994a), and some of the results from this work is included with changes
in notation.

Given a continuous stress history 0(t) with tmin S t S imam and the initial marten-
sitic fraction 5 = $0 at t = tmz-n, we need to find the resulting martensitic fraction
5 = g(t) within this time interval. For the purpose of modeling transformation
subloops, we divide the time interval into subintervals such that within each subin-
terval tk S t g t k +1, the stress is either nondecreasing or nonincreasing. Hence a list
of time sequences can be obtained as tmm = t1 < t2 < -- ' < t K :2 tmax. Although
such a division may not be unique, the kinetic rules used guarantee that the actual
representation of the transformation behavior is the same. Therefore, within any time

interval [tk, tk+1], we have either ('7 _>_ O or (’7 S 0.

22

If a 2 0 in time interval [tk, tk+1], the material undergoes forward transformation

and

_ _ 1-g(tk) _ FTU
g(t)—1 l—gFT(o(tk),T)[1 g ((t),T)]- (2-6)

 

Note here we only consider an isothermal situation, hence T does not change during
the analysis, it only affects the threshold stress for phase transformation. Similarly, if

('7 _<_ O in time interval [tk, tk+1], the material undergoes reverse transformation and

 

_ g(tk) RT 0
g(t)—gRT(0(tk),T)9 ( (t),T)- (2-7)

The above expressions implicitly use the assumption that at a given temperature the

threshold stresses for initiating forward and reverse transformation and completing

forward and reverse transformation remain constant. This was shown in the exper-

iments by Tanaka et al. (1995a). At a given temperature, the threshold stresses for

the subloops are the same as those for the major loop: forward transformation always
FT FT

starts at as (T) and completes at a f (T), and reverse transformation always starts

at 0§T(T) and completes at 0?T(T).

23

Substituting equation (2.4) into equation (2.6) gives equations for forward transfor-

mation associated with subloops as

 

 

 

93 0 S 0§T(T)i
9FT(0,93,T) = < 1— 1- gs{1+ cos [a — 05T(T)]W } oFT(T) < o < 0FT(T)
2 office — 05TH) 3 f
1 o > ofT(T),
\

in which gs is defined by

90 00 S 05TH"),

U “OFT 7f —
as = l 1_ 2(1— go){1+ cos #190) 8_ aga} 05m“) < 00 < gyro"),

 

 

where go = g(tk) and 00 = 0(tk) are the martensitic fraction and stress at tk when
the forward transformation interval starts. Note there is no explicit time dependence
in the expression of fiFT function, and the history dependent behavior is reflected in
its argument gs. The forward transformation function for major loop can be retrieved
by setting gs = 0 in equation (2.8). Figure 2.3 shows a schematic plot of austenite to

martensite transformation from austenite/martensite mixtures.

24

 

 

0.54

 

 

 

0'

 

 

6%) 05%; cm) can

Figure 2.3: Partial forward transformation curves. Regardless of the existing marten-
sitic fraction partial forward transformation from austenite to martensite always

starts at USFT(T) and completes at offrT(T).

Similarly, substituting (2.5) into (2.7) gives equations for partial reverse transforma-

tion associated with transformation subloops as

 

0 a < offeT(T),
~RT [0 - URT(T)l7r
g 0,9 ,T) = 9_8 _ f RT RT
( s 2{1 COSU§T(T)—0?T(T)} of (T)<o<os (T),
93 0' > 0.57171),
(2.10)
in which g3 is defined by
0 00 < ofr(T),
[0'0 -- awn

 

93 = 290{1— cos

90 00 2 0§T(T)-

25

 

 

(l
\

is

0.5-

 

/ 0'

are) 65%) 0.5% 6157‘s)

 

A

Figure 2.4: Partial reverse transformation curves. Regardless of the existing marten-

sitic fraction, partial reverse transformation from martensite to austenite always starts

at 0§T(T) and completes at 0?T(T).

Again 90 = g(tk) and 00 2: 0(tk) are the martensitic fraction and stress at tk when
the reverse transformation interval starts, and there is no explicit time dependence
in QRT function. The major loop for reverse transformation can be obtained by
setting gs = 1 in equation (2.10). Figure 2.4 shows a schematic plot of martensite to
austenite transformation from austenite/martensite mixtures.

Figure 2.5 shows the transformation subloops under cyclic loading. When stress os-
cillates between two fixed values, it will have converging subloops as the load cycle
increases. Figure 2.6 shows the stabilization of subloops under cyclic loading his-
tory. For more detailed discussions on hysteresis loops, please refer to Ivshin and

Pence (1994a, 1994b) and the stabilization of pseudoelastic loops is reported in the

experimental investigation by Wolons et al. (1998).

26

 

0.5 ‘

 

 

 

 

 

./ a; i .
0}”(T) 657T) 657T) 6}”le

0'

 

Figure 2.5: Transformation subloops under cyclic loading.

1.?) _____________

 

 

0.5 _

 

 

 

 

 

O L/ 1 l G
0}"(T) chlT) 053T) Gilli”)

Figure 2.6: Stabilization of subloops under cyclic loading.

2.2 Multi-dimensional Stress-Strain Relation

For the multi-dimensional setting, an equivalent stress is used to fit the above frame-
work in transformation kinetics. Here Von Mises stress (J2) is used as the effective
stress that drives the martensitic transformation, and the martensitic volume fraction
is determined by history of effective stress. To determine the components of the trans-
formation strain, additional assumptions are required. The transformation strain is

assumed to be co—axial to the stress deviator, and the magnitude of the transformation

27

strain (use some sealer measures) is assumed to be proportional to martensitic volume
fraction. With these assumptions, the transformation strain is fully determined by
martensitic fraction, calculated by (2.8) or (2.10), and the direction of stress deviator.
The total strain is

l

s = £6 + a". (2.12)

The elastic strain is

 

el 1 /\
= — — t 1 2.1
e 2,uT 2n(3/\+2n) r’r ’ ( 3)

where 7' is the stress tensor, p and A are two Lame’s constants. The transformation

strain can be expressed as

”-303
6 -20

S, (2.14)
where o is effective stress defined by a = V3.3 : S / 2 and S is stress deviator defined
by S = 7' — tr1'1. The evolution of E is the same as that of the one-dimensional
setting except that the dependent variable stress is now changed from uniaxial tensile
stress to effective stress. From equation (2.14), tretr = 0, indicating that the phase

transformation by itself is volume preserving in the present small strain setting.

One convenient scalar measure of the transformation strain magnitude is the effective

e" = ‘/ get?" : a". (2.15)

Immediately with equation (2.14) we have 6" = aé, implying that the magnitude of

strain et'r defined by

the transformation strain defined in equation (2.14) is proportional to the marten-

sitic volume fraction. Under uniaxial load, 7' = oel 8) e1, according to (2.14), the

28

resulting transformation strain is 8" = a5(e1 (3) e1 — 1/2e2 <8) e2 — 1/2e3 (8) e3). In e1
direction, it is easy to verify that the total strain coincides with equation (2.1). The
expression (2.14) is a deformation version of the flow rule in the model of Bondaryev
and Wayman (1988). Namely, their model gives the specification of detr rather than
a". The present deformation version follows closely the recent formulation of Briggs
and Ostowski (2002a) . It is expected that J2-style theories would generally involve
either a formulation like (2.14) or else an associated rate version (see for example
(3.25) of Qidwai and Lagoudas (2000)). Recently, for example, Briggs and Ponte
Castaneda (2002) have employed such a model to estimate the effective behavior of
shape memory alloy fiber composites.

The effective stress can be regarded as a scalar representation of the shear character
resident in the stress tensor. That this alone dictates the martensite phase fraction is
consistent with experimental results showing that shear is the predominant factor in
stress-induced martensitic transformation (Delaey et al., 1974). The above definition
of 5 and transformation strain imply tension / compression symmetry in the model un-
der consideration here. Although such models are standard, shape memory materials
generally display some tension/compression asymmetry (Gall et al., 1998 and Lim
and McDowell, 1999). A possible refinement of the model is to allow 5 = g(o, J3),
where J3 = detS as discussed for example in Gillet et al. (1995).

The advantage of J2 deformation theory rests on its simplicity. Although in com-
mon practice, J2 deformation theory is used only for proportional loading, Budian-
sky (1959) suggested that J2 deformation theory may be applied to nonproportional

loading paths under certain circumstances. In the context of plasticity, he argued that

29

the stress and strain in deformation theory also satisfies Drucker’s requirements for
the loading paths moderately deviated from the proportional paths for some yield sur-
faces. Therefore, in our study, the constitutive model modified from J2 deformation
theory is used to analyze both proportional and nonproportional loading paths.

Using the reverse transformation function defined by equation (2.10) and (2.11), and
the transformation strain defined in (2.14), the unloading behavior can be modeled.
An important distinction between this pseudoelastic model and the conventional plas-
ticity theory is that in conventional plasticity the plastic strain does not change when
the material undergoes unloading, while in the above model, regardless of loading and
unloading, the transformation strain (counterpart of the plastic strain) can change

due to the change of martensitic volume fraction, direction of stress deviator, or both.

30

Chapter 3

Plane Stress Analysis of Shape
Memory Annular Plates Globally

Loaded from A Phase

Starting from this chapter, we will focus on finding the responses of shape memory
annular plates subject to traction loads. The annular plate under study has an inner
radius ri, an outer radius r0, and a constant thickness (Figure 3.1). It is held at a
constant temperature above Af so that in the absence of mechanical loads the plate
is everywhere in the austenite phase. It is subject to uniform normal tractions p, on
its inner boundary and p0 on its outer boundary.

The tractions Pi and p0 are given by two continuous functions defined in time interval

tmz‘n _<. t S tmaa: 33
m = 15W),
P0 = 130(15)-

31

 

Figure 3.1: Generic phase distribution within the shape memory alloy annulus. The
inner ring M (Ti 3 r g rM) consists of martensite. The outer ring A (72A 3 r 3
r0) consists of austenite. The intermediate ring (‘3 (TM < r < TA) consists of an
austenite/martensite mixture.

We stipulate 1300mm) = Raltmin) = 0 so that at t = tmin we have an austenite
plate at t = th-n. Also, we assume that the traction loads change slowly such that
the effect of inertia on the deformation and phase transformation can be neglected,
therefore we can perform a quasi-static analysis. Although the responses of the plate
under the traction history are functions of time, the time scale does not affect the
final responses of the plate.

For a general traction history, the martensitic evolution at any point can follow either
the major loop or subloops as described in Chapter 2. The stress-strain relation is
load history dependent, which is not favorable in solving the boundary value problem
with an analytical approach. The history-dependency in the stress-strain relation
can be avoided given that the plate undergoes forward transformation starting from

austenite and there is no reverse transformation activated. In this case, starting

32

from an austenite plate, as time proceeds, the effective stress is non-decreasing ev-
erywhere in the plate. We use a global loading A plate, or a plate globally loaded from
A phase, to represent such a scenario. Another scenario includes the plate under-
goes reverse transformation starting from martensite without activation of forward
transformation. In this case, starting from a martensite plate, as time proceeds, the
effective stress is non-increasing everywhere in the plate . We use a global unloading
M plate, or a plate globally unloaded from M phase, to represent this scenario. We
use “global” or “globally” to emphasize loading (or unloading) throughout the plate.
Also, two qualifying phrases, “from A phase” and “from M phase”, are necessary for
the transformation to follow the transformation major loop.

Hence we begin by studying the global loading of A plates. This limits the analyses
to be valid for only a subclass of traction histories. In this chapter, the plane stress
analysis of shape memory plates globally loaded from A phase is performed; in Chap-
ter 4, the plane strain analysis of shape memory plates globally loaded from A phase
is performed. Finally, in Chapter 5, we will study the plane stress responses of shape

memory annular plates subjected to arbitrary traction histories.

3. 1 Problem Statement

Consider a thin plate subject to uniform normal tractions, the boundary conditions

of which are given by

Urr = —P0 at 7' = To,

(3.2)

arr = —p,- at r = Ti-

33

Note here we are performing a quasi-static analysis, therefore, the time parameter is
not important and it is omitted for conciseness. There is no shear traction on the
boundaries. This setting is consistent with axisymmetric radial displacement 217(7)
and vanishing azimuthal displacement Hg = 0. A plane stress solution of the plate
involves only arr and 099, both functions of r only.
The equilibrium equations reduce to the single equation

do" + a" _ “99 = 0. (3.3)

dr r

Attention is restricted to infinitesimal strains. The polar coordinate system provides

a principal frame in which the in—plane strains are given by

d'U/r Ur
Err —— F and 560 —— T. (3.4)

The out-of-plane strain 522 need not enter the problem formulation, but can be
determined from the constitutive description as shown in what follows. To maintain

compatibility, the in—plane strain components in (3.4) must satisfy

 

d e

34

The plate is composed of a shape memory material, and from equation (2.12) to

(2.14), the strains in terms of stresses can be written as

U '— VO'
Err = __7’L_ET__9_6 + 557‘,
— VO’

5‘ = —l/(O'rr + 006)
L zz E

 

 

+ at;

where E is the Young’s modulus, V the Poisson’s ratio, and eff, 53‘ and a?“ the trans-
formation strain components. Conventional thermal expansion is not considered here
as the analysis is isothermal. The transformation strain is obtained by specializing

(2.14) to the present case of plane stress as

tr _ 201T " 066

 

5r — 20 06,
55" = 20092; Orr 06, (3,7)
Etzr : _0902:0rra§
where o is the effective stress given by
a = 0,2,. — Orrdog + 036. (3.8)

As for this chapter, we restrict attention only to the forward transformation starting
from pure austenite and the case of isothermal loading at T > A f. For any given

traction load Pi and p0, we assume that there is a corresponding traction history

35

under which the plate is globally loaded from austenite to this traction load and the

reverse transformation is never activated. The existence of such a load history will

be discussed at the end of this chapter. For simplicity, temperature can be dropped

from the notation and only the forward transformation is considered. Therefore, only

(2.4) is involved in our analysis and it reduces to

 

0 fora g as,
1 _
9(0) = —[1 — cos( a 0‘3 7r)] for as < a < (If,
2 Uf—Us
1 foro 2 0f

by employing the following replacement:

gFT<a.T>—»g<a>. afT<T>—ws. afT<T>—+of.

The martensitic volume fraction is therefore uniquely determined by a as:

6 = 9(0)-

Immediately from (3.9), 5 = 0 if 0 S 03 and 5 = 1 if a 2 of.

(3.9)

(3.10)

(3.11)

In the present development, the martensitic distribution, 5 (r), is a continuous function

of r obeying 0 S 5 (r) g 1. For fixed p,- and p0 the plate can be in one of three broad

phase states: a fully austenite state denoted by A meaning that 5 (r) is identically

zero, a fully martensite state denoted by M meaning that 5 (r) is identically 1, or

some kind of combined state denoted through the use of 6 meaning that there is at

36

least one value of r such that 0 < 5 (r) < 1. The combined states can in turn be
organized into four types as follows: type (3 meaning that 0 < 5 (r) < 1 for all r; type
AC meaning that pure austenite is present so that 5 (r) < 1 for all r with at least one
value r such that 5 (r) = 0 and at least one other value r such that 5 (r) > 0; type
CM meaning that pure martensite is present so that 5 (r) > 0 for all r with at least
one value r such that 5 (r) = 1 and at least one other value r such that 5 (r) < 1; type
ACM meaning that both pure austenite and pure martensite are present such that
5(r) = 0 for at least one value 7‘ and 5 (r) = 1 for at least one other value r. The
possible phase distributions A, M, (‘3, AC, CM, ABM will be correlated to load values

Pi and pa in the ensuing development.

3.2 Reformulation of Governing Equations

Equations (3.6), (3.7) and (3.11) lead to the following expression for the in-plane

 

 

strains
8 —1(0 —1/o)+[ 1 —1l(a —-1-o)
rr - E rr 00 E3(0) El rr 2 06 1
(3.12)
5 -1(0 -V0‘ )+[ 1 —1](o —lo)
00 - E 09 77' ES(0) E 66 2 7'7.
where E3(o) is the secant modulus given by
E (0) — 1 (3 13)
S —1/E+ozg(o)/o' '

The first terms in the right hand side of (3.12) are the elastic strain, while the

second terms arise from the martensitic transformation. Substituting (3.12) into the

37

compatibility condition (3.5) gives

2 _.
060 0.71 138(0) = 0, (3.14)

1(a +0 )—
dr Tr 6“ 2E3(o) dr

in which the equilibrium equation (3.3) has been used to eliminate the terms involving
the Poisson’s ratio. The absence of the Poisson’s ratio is not unexpected for the case
of plane stress as discussed for example in Budiansky (1958) in the context of what
is known as the extended Michell’s theorem. This theorem gives that the stress field
for a traction boundary value problem is independent of Poisson’s ratio for the type
of problems under consideration here.

From equation (3.13) we have

 

d aE2(g — og’)
dr 3( ’ (Eag + a)2 ( ’
From equation (3.8) we have
do = 20” - 099 dom- + 2090 — arr dogg. (3.16)

 

 

d7 20 dr 2o dr

Substitute (3.8) and (3.16) into (3.14), use equation of equilibrium (3.3), dogg/dr can

be solved. And we obtain the following standard form of a two variable ODE system:

 

dUrr 0rr - 009
dr 7' ’

(3.17)
(1099 _ 3(0rr1009)0rr - 009

dr A(orr, 099) r

 

1

38

where

A(0r,~, 099) = 403 + Ea[30,?rg + og'(0rr — 2099)2],
(3.18)

B(07~r, 099) = 403 + Eoz[402g + (g — ag')(202 — 30rr099)l,

and for simplicity g(0) is written as g and g' = g’ (0) denotes the derivative dg/d0.
Note that A(0,~,~, 099) > 0 as long as 0 7é O, and this equation can be use to integrate
for a solution if we know arr and 099 at any given r value.

From equation (3.9), the effective stress dictates the evolution of martensitic fraction.
However, it not obvious from (3.17) how the effective stress distributes, therefore it
is desirable to reformulate equation (3.9) in terms of effective stress. This is accom-

plished by introducing the following change of variables due to Nadai (1950)

2
Urr = — sin(fi +
fl (0 > O) (3.19)

I 2
l 009 = 75 sinw — 3,90.

where 6 is a function of r and for convenience is taken as obeying 0 S E < 27r. This
transformation is valid when 0 79 0. When 0 = 0, both arr and 099 are necessarily
zero. In terms of 0 = 0(7) and [3 = E (r), the equilibrium (3.3) and the compatibility

(3.14) become

 

 

sin(fl+ (3)3— +ocos(B+ 6)%—5+ acisfi = 0,
(0 > 0) (3.20)
sinfl— +0cos fid_fi _ asin(fi— 7r/3)dE3(0) = 0.

2E3 (0) d?"

39

By virtue of the definition of Es(0) given in (3.12) the above equations can be written

as

 

 

Sinw + 25%:- + acosm +9515? ._. _Ucisfi,
(0 > 0) (3.21)
. Ea(9 - 09’) . 7r d0 dfi _
[smfi _ 2(Eag +0) 3mm — 3)]3‘; + 000353; _ 0,

Solving the system (3.21) for d0/dr and dfi/dr gives

d0 __ 02 0052 B

 

 

 

E: _ T’
(0 > 0) (3.22)
— = isni — €321.53? ssii - a]
with A given by
A = %[1 — Egggz—Uag') cos2(fi + g)]. (3.23)

It is useful to note that the dependence of governing equations (3.17), (3.18), (3.22)
and (3.23) on parameters E and a is solely through the product Ea.

By definition 0 > 0 and E, a, g and g’ are non-negative, therefore A > 0. It follows
from (3.22) with A > 0 that

_ < 0. (3.24)

Thus the maximum effective stress amax occurs at r = r,- and the minimum effective
stress 0min occurs at r = r0. Accordingly 5 (r) 18 also a non-increasmg functlon of r,

i. e. d5 / dr 3 0. The six possible phase distributions: A, M, (3, AC, BM, ACM, now

40

correspond to

 

A: 5(ri) = 0 <=> 0(7‘2) S 03,
AB: 1> 5(r)) > 0, 5(r0) = 0 41) of > 0(ri) > 03, 0(r0) S 03,
ABM: 5(rzi) = 1, 5(r0) = 0 4:) 0(r,-) 2 0f, 0(r0) g 03,
< (3.25)
G: 1> {(73) > O, 5(r0) > 0 4: of > 0(r,-) 2 0(r0) > 03,
BM: 5(r,-) =1, 1> 5(r0) > 0 4i) 0(ri) 2 0f, of > 0(r0) > 03,
( M: 5(r0) =14: 0(r0) _>_ 0f.

The phase distribution AB involves an austenite outer ring on r A < r < r0 and a
mixture zone on an inner ring r, < r < r A The separating interface 7‘ = r A is
associated with satisfaction of the condition 0(rA) = 03. The phase distribution BM
involves a mixture zone on an outer ring TM < r < r0 and a martensite inner ring
rz- < r < rM. The separating interface r = TM is associated with satisfaction of
the condition 0(rM) = 0 f. The phase distribution ABM involves three rings: an
austenite outer ring 7' A < r < To, a mixture zone on an intermediate ring TM < r <
r A, and a martensite inner ring r, < r < TM: These are the same type of structures
as obtained by Birman (1999).

From the first equation of (3.22), it follows that a constant effective stress field requires
either 5 = 7r / 2 or B = 37r / 2. These correspond to uniform equi-biaxial stress solutions
0rr(r) = 099(r) = —p so that 0 = Ipl. Such solutions occur if and only if p0 = 171' and
are the only solutions with effective stress 0(r) = const. Here p,- = p0 = p > 0 gives

B = 37r/2 and p,- = p0 = p < 0 gives 6 = 7r/2. The phase distribution for uniform

41

equi—biaxial stress is respectively of type A, B, M according to whether |p| g 03,
on < Ipl < 0). Ipl 2 of.

For the generic case p0 7S p,, it follows that (3.24) holds with strict inequality. To
obtain the effective stress distribution requires integration of (3.22). Note for example
that, if {0*(r), fi*(r)} is a solution of (3.22) corresponding to boundary conditions
Orr = —p;‘ at r = rz- and arr = -p(*, at r = r0, then {0*(r),B*(r) + 7r} is a solution
corresponding to boundary conditions O'rr = +1): at r = r,- and 077- = +19; at
r = r0. Therefore in the (pi, p0) plane, two load pairs that are related to each other
by a 180° rotation about the origin will produce the same distribution of effective
stress and hence the same phase configuration in the plate. From equation (3.17),
we can replace the Orr and 099 with —0rr and ~099 respectively while maintain
the equality of the equation. Such equality reiterates the symmetry property in the
(pi, p0) plane. This symmetric property is inherent in the constitutive model in view
of its tension / compression symmetry.

If arr and 099 are specified at some location r, then equations (3.17) can in principle
be integrated for arr = 077-(r), 099 = 099(r). However, boundary conditions in the
original form (3.2) do not generate conditions for (Orr, 099) at any one point. Gener-
ally, the shooting method, as described for example by Roberts and Shipman (1972),
can be used in the case of two-point boundary conditions such as (3.2). Starting
from 099(ri) given by the the elastic solution for p,- and p0, we can perform a root-
searching, such as Newton’s iteration, to find the correct 066(7‘2') so that integration

of equation (3.14) results in 0T7-(r0) = —po.

42

Similarly, if 0 and B are specified at some location 7‘, then equations (3.22) can in
principle be integrated for 0 = 0(r), B = 3(r). Again the shooting method can be
used to obtain the solution associated with the two-point boundary conditions (3.2).
The shooting method can be avoided for the case of alternative boundary conditions
wherein both 0 and arr are specified at one boundary, in which case (3.19) is used to
calculate the boundary value [3 (modulo an inconsequential trigonometric periodicity).
The specification of both 0 and arr at a common boundary is motivated by (3.25)
for the special case when the specified value of 0 is either 0 f or 03. In particular this
provides a direct procedure for studying the four special conditions r A = ri, r A = r0,
TM = ri, TM = To that are associated with transitions between the different phase

distribution types in (3.25).

3.3 A Fully Austenitic Plate (A Phase Distribu-
tion)

This section considers plates that are completely in the austenite phase and so corre-
spond to a phase distribution of type A in (3.25). This requires 0 S 03 everywhere

in the plate, hence g = g’ = 0, A = 0 / 2 and equation (3.22) reduces to

do 20coszfi d5 2sin [3 cos 6
——-—-——— and -——=——————.

_ _ .2
dr 7' dr 7‘ (3 6)

43

The solution to equation (3.26) is given by

osinfi = cl and sin fi/ cosfi = c2r2 for cos ,6 51$ 0,
(3.27)

0 =03 for cosfi=0
where cl, c2 and C3 are constants to be determined by the boundary conditions.
Radial stress arr and hoop stress 099 can be obtained by substituting (3.27) into
(3.19) as

for cosfi 74 0,

Cl
0'7”]. —— C]. + _—

Cl and 0 - c
\/§C2T2 00 1
ans = 099 = :l:C3 for cos B = 0

and it follows from this equation that cos B = 0 if and only if p,- = p0.

If p,- # p0, then the constants c1 and 02 can be obtained from the first equation of

 

 

(3.28) as
2 2 2
p 7‘ par p'r- -por
z 22 2 0 and c2 = z z 02 2 for p,- #po. (3.29)
7'0 - Ti £070 _ Pilri To

If instead p,- = p0 = p, then the second equation of (3.28) gives C3 = |p| whereupon

0rd?) = (7990‘) = —p, 0(7‘) = Ipl for Pi = pa. (3-30)

44

Equations (3.28), (3.29) and (3.30) together retrieve the classical linear elastic solu-

tion, with in-plane stresses given by

 

 

 

 

 

 

a _ W3 - 19ng (p0 - pad-4%
7‘7“ — 1
r3 — 7'? (n2, — T32)T2 3 31)
. 2 _ 2 _ . 2 2 ( °
096 _ 1927} para __ (P0 Pilri To
rg - r22 (r3 — 7‘z.2)r2
The effective stress is given by
A \/(7‘c22100 — ifs»? + Sigism- — po>2/r4
0 = 0(r) := 2 2 . (3.32)

’"0 ‘0'

The maximum and minimum effective stresses are given by Umax = 0(ri) and 0min =
0(r0) respectively.
The validity of this solution requires Umax S 03. The case amax = 03 is associated

with parameter values obeying

(r;1 + 37:3)pz2 — 27:30"? + 3rg)pipo + 47/3193 = (r3 — r22)203. (3.33)

Equation (3.33) defines an ellipse in the (pi, p0) plane containing the origin (Pi: p0) =
(0,0). Values (pi, p0) internal to this ellipse are associated with an A phase distribu-
tion, while values (103', p0) outside this ellipse involve one of the other phase distribu—
tions in (3.25). This ellipse will be referred to as the Fully Austenitic Loop (FAL) in
the (pi, p0) plane and corresponds to the transitional condition r A = r,. The FAL

is symmetric with respect to 180° rotation about the origin, a consequence of the

45

tension/ compression symmetry in the modeling description. In the sequel additional
curves will be constructed in the (pi, p0) plane, thus giving rise to a fully articulated
structure map that distinguishes between the phase distribution possibilities (3.25).

An alternative way of obtaining the PAL is to solve (3.26) subject to a requirement
”(Til = 03. If 6(a) = fig 75 7r/2,37r/2, then according to (3.27) CI and 02 are given

by c1 = as sin [30 and c2 = tan fio/rzz, hence the stresses are given by

or,» = as [sin 60 + C———OS\/§O;: ],
099 = as [sin fio— 13:8;ng ], (3.34)

 

4
7'.
0 = 03 \/sin2 50 + cos2 50 —Z.
r

Note this solution is also valid for 60 = 7r/ 2, 3r /2. Applied tractions p,- and p0 are

then given by

p,- = —(sin 30 + (1ng )03 and p0 = —(sin 00+ +C————(:S/§f;zz. )0 (3.35)

 

Equation (3.35) defines an ellipse in parametric form within the (pi, p0) plane where

60 is the varying parameter. Solving (3.35) gives sin [30 and cos 50 as

 

 

2 2
-r- — r
sin 60 2 p22, 12% 0 and cos BO: [(190 —2p,)r0’ (3.36)
(’0 “ 7‘; )Us (To " 7'2: )03

46

whereupon the identity sin2 fig + cos2 60 = 1 retrieves (3.33). The use of 6 as a
varying parameter in obtaining loops associated with the remaining transitional cases

7‘ A = r0, TM = ri, TM = To is central to the remaining development of this Chapter.

3.4 A Fully Martensitic Plate (M Phase Distribu-
tion)

This section considers plates that are completely in the martensite phase and so
correspond to a phase distribution of type M in (3.25). Hence g = 1 and g’ = 0, and

equation (3.22) reduces to

 

 

dz _ _02 cos2 6
dr _ rAM ’ (3 37)
d6 _ acosfi . Ea . W
H? _- rAM [smfi 2(Ea + 0) sm(fl 3)]
where AM is obtained from (3.23) as
AM = Z [1 — Ea cos2(,6 + Z)]. (3.38)
2 Ea + 0 6

A closed form solution to (3.38) is not obvious, hence prompting a numerical approach.
Such solutions are valid provided 0(r0) Z 0 f.

For prescribed (171:, p0), the solution to (3.38) can be obtained by the shooting method.
As an example, Figure 3.2 displays the radial variation of 0, ans and 099 for a plate

with E = 50 GPa, a = 0.05, 03 = 75 MPa, 0f = 300 MPa, ro/rz- = 2, p0 = 503 and

47

 

 

 

 

 

 

 

Figure 3.2: Typical stress distributions (normalized by 03) for a plate with M phase
distribution (0(r) Z 0 f). Here E = 50 GPa, 0f = 403 = 300 MPa, ro/ri = 2 and
(pi, p0) = (1, 5)03. Solid lines represent the case of a = 0.05. For comparison, stresses
for an elastic solution ( = 0) under the same boundary conditions are plotted as
dashed lines.

p,- = 03. The corresponding stress distributions for an elastic material (a = 0) subject
to the same normal tractions p, and p0 are plotted as dashed lines. Note that the 0m
distributions for these two materials are similar, but the 099 distributions are quite
different. In particular the two materials give different ratios 0rr/099 therefore raising
questions with respect to any assumption that such a ratio is preserved (Birman,
1999).

The equality case for the validity condition 0(r0) 2 0 f defines the Fully Martensitic
Loop (FML) in the (pi, p0) plane. The FML at fixed 0 f, Ea and rO/ri can be obtained
from integrating equation (3.37) with boundary condition (0(r0), fi(r0)) = (0 f, 60)
upon varying 60 from 0 to 27r. Like the FAL, the FML is also symmetric with respect

to 1800 rotation about the origin. The FML corresponds to the transitional condition

48

TM = r0, and the region exterior to the FML corresponds to tractions (pi, p0) that
generate an M phase distribution. The FML will surround the FAL in the (pi, p0)
plane, although unlike the FAL, the FML is generally not an ellipse. The region
between the FAL and the FML corresponds to loads (pi, p0) that generate one of the
phase distributions containing a mixture zone: B, AB, BM or ABM.

Although the F ML is usually determined numerically, it will pass through (pi, p0) =
:i:(0 f, 0 f) which provides the uniform equi-biaxial stress solutions on the FML. More
generally, the line p,- = p0 provides the uniform equi-biaxial stress solutions in the
(pi, p0) plane. On this line 0rr(r) = 099(r) = —p, giving a uniform effective stress
0(r) = |p| and hence a phase fraction field 5 (7) that is also independent of r and
determined directly from p as 5 = g(lpl). Hence such loading results in a phase
distribution of either type A, B or M.

Figure 3.3 displays both the FAL and the FML for the parameters E = 50 GPa,
a = 0.05, 03 = 75 MPa, 0 f = 300 MPa, and ro/rz- = 2. Note that the region enclosed
by the FML is not convex for these parameters. Note also that the FML is highly
elongated in the pi direction compared with the po direction. In this case the external
traction p0 is significantly more efficient than the internal traction pa in promoting
an M phase distribution. For example, a traction free inner boundary requires an
outer traction p0 = 3.3903 to fully martensitize the plate, whereas a traction free
outer boundary requires an inner traction p,- = 27.7403.

An FML is completely determined by the parameters 0 f, Ea and ro/ri. Figure 3.4
shows the effect of varying 0 f on the F ML at fixed Ea and ro/ri. In particular, this
figure displays a family of seven FMLS corresponding to 0 f /03 = 1,2,3,4,6,8,10,

49

Pb/os

8

6

4

2’ .

. pi/U
/4 . 1 . 1M 1 1 Is

 

 

L l
40 —3o’ -20 -10 0. 10 20 3o 40
v-
-6 L
-3 '_
-10 -

Figure 3.3: The FML (Fully Martensitic Loop) and FAL (Fully Austenitic Loop)
for E = 50 GPa, a = 0.05, 0 f = 403 = 300 MPa and ro/rz- = 2. Tfactions inside
the smaller ellipse (FAL) produce a fully austenitic plate, while tractions outside the
bigger loop (F ML) produce a fully martensitic plate.

with the other parameters set at E = 50 GPa, a = 0.05, 03 = 75 MPa and ro/ri =
2. As one would anticipate, the F MLs are nested with respect to increasing 0 f.
The small internal ellipse is the common FAL. Note that the FML corresponding to
0 f / 03 = 1 intersects the FAL at the uniform equi-biaxial stress solutions (pi, p0) =
i(0 f, 0 f). Thus if 0 f = 03 then equi-biaxial loading such that p,- = p0 = p results
in a uniform phase distribution of type A transitioning directly to a uniform phase
distribution of type M as the common traction value p passes through the common
threshold stress 0 f = 03.

Figure 3.5 shows the effect of varying Ea on the FML at fixed 0 f and ro/ri. In partic-
ular, this figure displays a family of eight FMLS corresponding to Ea / 530 MPa = 0,
0.5, 1, 2, 3, 5, 7, 10, with the other parameters set at of = 403 = 300 MPa, and

ro/rz- = 2. The common FAL is again shown as the small internal ellipse. All

50

po/O's
15¢

lO/\10

fiz pi/O's

. /: I . 1 . —_ .
-50 40 w 10 5 20 1 40 50

-10

 

 

-15

Figure 3.4: FMLS for various 0 f values. Here E = 50 GPa, a = 0.05, 03 = 75 MPa

and ro/rz- = 2. The smallest loop is the common FAL, and the other curves are
FMLS, each marked with corresponding value of 0 f / 03.

the eight FMLs pass through the two common uniform equi-biaxial load points
(pi,p0) = :l:(0f,0f). Otherwise the FMLS are separate and elongate significantly
along the p,- direction with increasing Ea. For Ea = 0 equation (3.37) reduces
to (3.26), whereupon the classical elastic solution can be used to plot the FML for
Ea = 0 as an ellipse in the (1),, p0) plane.

Figure 3.6 shows the effect of varying ro/rz- on the FML at fixed 0 f and Ea. In
particular, this figure displays a family of seven FMLS corresponding to ro/ri = 1.1,
1.3, 0.5, 2.0, 2.4, 3.0, 10.0, with the other parameters set at E = 50 GPa, a = 0.05
and 0 f = 403 = 300 MPa. All the seven FMLS pass through the two common
uniform equi-biaxial load points (pi, p0) = :l:(0 f, 0 f). Otherwise the FMLS elongate
significantly along the p,- direction with increasing ro/ri. There is also a much less

pronounced contraction of the FMLS in the po direction with increasing ro/r, so that

51

 
 
 

 

1’ 4 1 . 1 n 4 .

-50 -40 -30’ -20 110/

\ g /
low” 0 -

 

Figure 3.5: FMLs for various Ea values. Here 0 f = 403 = 300 MPa and ro/rz- = 2.
The smallest loop is the common FAL, and the other curves are FMLS, each marked
with the value of Ea/ (500 MPa).

intersection of any two such FMLs now also takes place off of the line of uniform equi-
biaxial stress solutions (pi, p0) = (p, p). In the limit ro/r, -—> 00 the FML asymptotes

to the two horizontal parallel lines p0 = :l:0’f.

3.5 A Plate Containing Mixtures of Martensite
and Austenite (Phase Distributions of type:

AB, ABM, CM and c)

We now turn to consider the remaining types of phase distribution given in (3.25),
namely types AB, ABM, BM, and B. The structure map is useful in this discussion, as
all such phase distributions occur in the region between the FAL and the FML. Here
it is useful to recall that the phase distribution in this region is ensured to be of type
B on the line of uniform equi-biaxial load (pi,p0) = (p,p), 03 < Ipl < 0 f. However

off of this line it is not yet clear what type of phase distribution occurs. To make this

52

 

 

 

ILL P
1
[fr .4\
-50 «'1 - 0 - 0 1 t 0 pi /0's
1

 

 

 

Figure 3.6: FMLS for various rO/ri values. Here E = 50 GPa, a = 0.05 and 0 f =
403 = 300 MPa. The smallest loop is the common FAL, and the other curves are
FMLS marked with ro/ri values.

determination it is useful to complete the structure map via the construction of two
additional curves.

The first such curve is the Austenite Disappearing Loop (ADL) that is associated
with solutions obeying 0(r0) = 03. Equivalently, the ADL is associated with the
transitional condition r A = r0. Austenite (unmixed with martensite) can only be
present at locations on the structure map inside of the ADL. Thus inside the ADL
the phase distribution must be of type A, AB, or ABM. Outside of the ADL the
phase distribution must be of type B, M, or BM. The ADL is strictly within the
FML so long as 0 f > 03. The ADL is also external to the FAL except at the two
points of intersection (pi, p0) = d:(0s, 03) which are on the equi-biaxial load line of
the structure map.

The final curve for completing the structure map is the Martensite Emerging Loop
(MEL) that is associated with solutions obeying 0(r,-) = 0 f. Equivalently, the MEL

is associated with the transitional condition r, = rM. Martensite (unmixed with

53

austenite) can only be present at locations on the structure map outside of the MEL.
Thus outside the MEL the phase distribution must be of type M, BM, or ABM.
Inside of the MEL the phase distribution must be of type B, AB, or A. The MEL
is strictly outside the FAL so long as 0 f > 03. The MEL is also internal to the
F ML except at the two points of intersection (pi, p0) = :i:(0 f, 0 f) which are on the
equi-biaxial load line of the structure map.

With the exception of locations on the equi-biaxial line, both the ADL and the MEL
are determined numerically. This is most easily accomplished using the same pro—
cedure that generates the FML upon replacing the FML condition 0(r0) :2 0 f with
either the ADL condition 0(r0) = 03 or the MEL condition 0(r,-) = 0 f. Like the
FML, both the ADL and the MEL are generally not elliptical. Like the FAL and
the FML, both the ADL and the MEL are symmetric with respect to 180° rotation
about the origin, hence the symmetry of structure map follows.

Two examples of fully articulated structure maps are presented in Figure 3.7 and
Figure 3.8. Both correspond to material parameters E = 50 GPa, a = 0.05, 03 =
75 MPa, 0 f = 300 MPa. Their difference is due to differing values of ro/ri, namely
ro/rz- = 2 for Figure 3.7 and ro/rz- = 5 for Figure 3.8. Each structure map is
partitioned into various regions by the four loops FAL, FML, ADL and MEL. Each
region corresponds to a distinct phase distribution type as also displayed in Figure
3.7 and Figure 3.8.

The qualitative difference between these two structure maps is that in Figure 3.7
the ADL is strictly inside the MEL while in Figure 3.8 it is not. We shall refer

to the former as a type I structure map and to the latter as a type 11 structure

54

 

 

 

 

 

Figure 3.7: Structure map (type I) for ro/rz- = 2. Here E = 50 GPa, a = 0.05 and
0 f = 403 = 300 MPa. Note that there is no ABM region. Traction pair (p,,p0) =
(5, 0.5)03, marked as a dot, produces a BM phase distribution. The resulting stresses
are shown in Figure 3.10.

map. A type II structure map has two symmetric regions corresponding to an ABM
phase distribution whereas a type I structure map does not involve this type of phase
distribution. To within unit scaling, the structure map is determined by the material
parameter combinations Ea, 0 f, 03 and geometry parameter ratio ro/ri. For fixed
material parameters Ea, 0 f, 05 one finds that the structure map is of type I for
small ro/rz- but is of type II for large rO/ri. The condition r0 sufficiently greater than
r,- makes possible a situation supporting the condition r, < TM < r A < r0 that is
associated with the phase distribution ABM.

For given Ea and 0 f, 03 there is a special value of ro/ro such that the structure
map is of type I below the special value but of type II above the special value.
For the material parameters associated with Figures 3.7 and 3.8 (Ea = 2.5GPa,
0 f = 403 = 300 MPa), this special transitional value is found to be ro/r, = 3.58.
The structure map for Ea = 2.5GPa, 0 f = 403 = 300 MPa, and ro/rz- = 3.58 is

plotted in Figure 3.9 showing how the ADL osculates the MEL.

55

 

 

 

 

 

 

 

 

 

pi/O's

Figure 3.8: Structure map (type II) for ro/rz- = 5. Here E = 50 GPa, a = 0.05
and 0 f = 403 = 300 MPa. All six types of phase distributions are present. Traction
pair (19,3190) = (5, 0.5)03, marked as a dot, produces a ABM phase distribution. The
resulting stresses are shown in Figure 3.11.

The stress distribution for the phase distributions under discussion can be determined
numerically on the basis of (3.22) using the shooting method, in a similar fashion to
this determination for a fully martensitic plate. For example, Figure 3.10 displays
the radial variation of 0, 0”, 099 for the BM phase distribution associated with the
load (pi,p0) = (4, 0.5)03 represented by the dot in Figure 3.78 (Ea = 2.5 GPa, 0f 2
403 = 300 MPa, ro/rz- = 2). A mixture of austenite and martensite appears on the
outer plate boundary r = r0. This mixture becomes progressively richer in martensite
as r decreases to r = TM = 1.39r,- within which the plate is in the martensite phase.
This solution could, in principle, be obtained by many quasi-static load paths in
the (pi,p0) plane so long as the effective stress at each point in the plate remains
non-decreasing. In particular, the proportional loading path (pi,p0) = (5,9.5)03k
where k is a load path parameter is associated with an initial phase distribution of
type A. There will then be three values of k < 1, say k1 < k2 < k3, associated

respectively with transitions through the sequence of phase distributions: A, AB, B,

56

 

 

 

BM

po/Gso

T“ I

-2 ' BM

 

 

 

 

 

 

-10 A 5 A (i A 5 [0

pi /0's
Figure 3.9: Structure map for the special value ro/rz- = 3.58. Lower values produce
type I structure maps, while higher values produces type II. Here E = 50 GPa,

a = 0.05 and 0 f = 403 = 300 MPa. There is no ABM region, and ADL contacts

MEL at two points. One of the two points (shown as a dot) represents (pi,p0) =
(3.565, 0.995)03, producing a B phase distribution with 0(r,-) = 0 f and 0(r0) = 03.
The resulting stresses are shown in Figure 3.12.

BM. Here k = k1 is associated with the emergence of the AB front r A at the inner
plate boundary. This front proceeds through the plate as It increases until it merges
with the outer boundary when k = [92. The BM front TM then emerges at the inner
plate boundary when k = k3 and continued loading drives this front to r = 1.24rz-
when k: = 1. Further loading will then drive this front to the outer boundary, resulting
in a phase distribution of type M.

Figure 3.11 shows the stress distributions at the same load (pi, p0) = (5, 0.5)03 as in
Figure 3.10 but for the structure map of Figure 3.8. In this case the larger value of
r0 (r0 = 57’3‘) gives rise to a phase distribution of type ABM with r A = 4.84r,~ and
TM = 1.23m. In this case the proportional loading path (Pi: p0) = (5, 0.5)03k passes
through phase distribution states of type A, AB, ABM with the transitions taking
place at two special k values that are analogous to k1 and k3 of the previous example

associated with Figure 3.10.

57

 

 

 

 

 

 

 

 

Figure 3.10: Stresses (normalized by 03) for a BM phase distribution. Here E =
50GPa, a = 0.05, of = 403 = 300MPa, ro/ri = 2 and (pi,po) = (5,0.5)03. An
elastic solution (a = 0) is shown for comparison (dashed lines).

A final set of stress distributions is shown in Figure 3.12. These correspond to the
osculation point of ADL with MEL at (p,, p0) = (3.565, 0.995)03. Along the propor-
tional loading line from the origin through this point, this solution is associated with
the disappearance of the pure austenite phase at the outer plate edge (r A = r0) that
is simultaneous with the appearance of the pure martensite phase at the inner plate
edge (TM = ri).

For the case of a type 11 structure map, Figure 3.13 summarizes the relationship
between the structure map and the effective stress distributions. Each marked point in
the structure map corresponds to an effective stress distribution graph. The points are
classified into two categories: those on the bounding loops, and those inside an area of
unique phase configuration type. The coordinates of the points are listed in Table 3.1.

The effective stress distributions are plotted in small scale figures from Figure 313(1)

58

 

 

m = 4.847,-

 

 

-3.0: _ __ 3.5 .. - .40.... . 5.0 T/Tz:

rr’ ””””
I

   
 

 

Figure 3.11: Stresses (normalized by 03) for a ABM phase distribution. Here E =
50 GPa, a = 0.05, 0f = 403 = 300 MPa, ro/rz- = 5 and (plypo) = (5,0.5)03. An
elastic solution (a = O) is shown for comparison (dashed lines).

to Figure 3.13(f) A similar summary can also be made for type I structure maps (as
in Figure 3.7) and for the non-generic structure map that transitions between type
I and type II (as in Figure 3.9). For the transitional structure map (as in Figure
3.11), points 3, 5, 6 and 7 merge to a common osculation of ADL with MEL, thus
annihilating the ABM region. For a type I structure map, points 3, 5, 6, 7 and c have

disappeared.

3.6 Solutions for Global Unloading of M Plates

For plane stress solutions of shape memory plates globally unloaded from M phase,
the evolution of martensitic transformation follows equation (2.5) everywhere in the

plate. Again the martensitic volume fraction is uniquely determined by the effective

59

 

 

 

 

 

 

 

 

.4-

Figure 3.12: Stresses (normalized by 03) for a special B phase distribution with
TM = r,- and TA = r0, for E = 50 GPa, a = 0.05, 0f = 403 = 300 MPa, ro/ri = 3.58
and (pi, p0) = (3.565, 0.995)03. An elastic solution (a = 0) is shown for comparison
(dashed lines).

stress 0. As a result, we can use the same expression (3.9) for martensite evolution

by setting 03 = 0;?!“ (T) and 0 f = 033T

(T). Since no change of other equations is
required, the solutions can be obtained in the same way as finding solutions for global

loading of A plates.

 

 

Point 11,- /03 po /03 Point p, /03 po /03
1 0.903 0.390 9 4.000 4.000
2 1.000 1.000 10 2.386 3.951
3 3.496 0.543 a —0.600 —0.200
4 2.287 1.055 b —2.000 —0.500
5 3.689 1.136 c -4.000 —0.600
6 3.206 —-.096 d ——3.000 -—2.000
7 5.226 1.154 8 —6.000 —2.000
8 4.294 2.605 f —6.000 —5.000

 

Table 3.1: The coordinates of the points on Figure 3.13.

60

The procedure of global unloading from M phase to A phase for a given shape memory
material can be viewed as the reverse of global loading from A phase to M phase for
a virtual material. For this virtual material, the threshold stresses initiating and
completing the forward transformation equal the threshold stresses initiating and
completing the reverse transformation for the given material respectively. Given time
t within interval [t0,t1], let {p,(t),p0(t)}, 5(r, t) and 0(r, t) be the corresponding
traction history, martensitic fraction and effective stress; then we have 5 (r, to) = 1
and 5(r,t1) = 0 as the plate is transformed from M phase to A phase. Global
unloading implies 0(r, ta) 2 0(r, tb) for to _<_ ta < tb 3 t1. For an A plate composed
of the virtual material subjected to load {At (t), p0(t)} = {15,- (to +t1 —t), 130(3) +t1 —t)}
(to g t 3 t1), the martensitic fraction and effective stress distribution will be 5 (r, t) =
5(r,t0 + t1 — t) and 0(r, t) = 0(r,t0 + t1 — t). Note that 0(r,ta) _<_ 0(r,tb) for
to g ta < tb _<_ t1, 5 (r, to) = 0, and 5(r,t1) = 1, therefore the plate composed of the
virtual material is globally loaded from A phase under traction {p‘z-(t), p0(t)}.

For a traction history that drives the plate globally unloaded from an M plate, we
can use unloading structure map to determine the phase combination. The structure
maps associated with the global loading of A plates will be referred to as loading
structure map. Similar to the loading structure map, four loops partitioning (pi, p0)
plane into regions with different phase combinations complete the unloading structure
map. These loops can be obtained using the same procedure as the construction of
the loading structure map by changing 05T(T) to off-(T), and 0§T(T) to a? (T).
Figure 3.14 shows the loading structure map and unloading structure map together in

(pi, p0) plane. The loading structure map is plotted as solid curves, and the unloading

61

structure map is plotted as dash curves. The material constants are: E = 50 GPa,

a = 0.05, of.” = 75MPa, agiT = 150MPa, of T = 225MPa, and of T = 300MPa.
Note that loading structure map is type II but unloading structure map is type I.
Consider the following traction history {152-(t), p0(t)}

13,-(t) = 2t00 and p0(t) = 4t00 for 0 S t _<_ 1,
(3.39)

13,-(t) = 2(2 — t)00 and 1300) = 4(2 — t)00 for 1 < t S 2,

where 00 = 0?T(T). Under this traction history, for an annulus with r0 = 27‘,- and

the above material, it is first globally loaded from A phase to M phase, and then
globally unloaded back to A phase.

The solid line with two arrows in Figure 3.14 is the trajectory of {15,-(t), p0(t)} defined
in (3.39). For 0 g t S 1, the annulus is globally loaded from A phase to M phase
because {132(1), 130(1)} 2 {200, 400) (point 2 in Figure 3.14) locates outside the FML
of the loading structure map. For 1 S t S 2, the annulus is globally unloaded from
M phase back to A phase. The reverse transformation will not begin until point
3 is reached. The global loading / unloading assumption is verified numerically by
checking the difference of the effective stresses for any given r between two consecutive
loading steps. Here 200 time steps with At = 0.02 and 501 points in radial direction

with Ar = 0.002 are used.

62

3. 7 Discussions

In this chapter, a boundary value problem is formulated and solved for the equilib-
rium stress and strain fields in an annulus composed of a shape memory material
which is globally loaded from A phase. The effective stress and the martensite phase
fraction are each shown to be a non-increasing function of r. This restricts the phase
partitioning to be one of six generic types (arranged along the direction of decreasing
r): A, B, M, AB, BM, and ABM. Other combinations violate either the continuity
of the martensitic fraction (such as AM), or the non-increasing of the martensitic
fraction (such as MBA).

The phase partitioning is described by a structure map in the (pi, p0) plane that pro-
vides a straightforward guide for finding loads that are required to produce a desired
phase partitioning. Parametric studies show how geometry and material properties
affect the structure maps. Such information would generally be useful for systems in-
volving device control by means of shape memory alloy plates and tubes. If warranted,
the strain and displacement fields are easily obtained from the stress distribution us-
ing equations (3.6) and (3.4).

Global loading from A phase, or global unloading from M phase, is the prerequisite to
determine the martensitic fraction simply by (3.9). Intuitively, a proportional loading
path should produce a global loading plate. Although this is also true in practice,
it is not obvious from the equations that such a property holds, and an analytical
proof is not available. Therefore, verifications of the global loading condition should

be performed. Nonproportional traction histories can also produce global loading

63

solutions. If a traction history produced a global loading solution, from continuity,
a neighborhood of the traction history should usually also produce a global loading
solution. However, under a general traction history, more often a plate is neither
globally loaded from A phase, nor globally unloaded from M phase. Further studies
on shape memory annular plates under arbitrary traction histories are developed in

Chapter 5.

64

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

J
4 - l
l
2 i- a
190/63 0 - -
-2 n i
l
-4 p _
l ' f
'6-8 -3 ‘ it -5 T 6 1 3: is 8
of of c’f of
Us 0'3 Us \ c’s ¥
0:, Cf Cf \\ Of
Us 03 C53 “3
Ti (5)8 7"0 Ti (6)8 7'0 7'7; (7) me 7'0 7'1; (8)M To
°f "IL “I “f

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

""1; (9)M To 7'3 (10)M 7'0 7'2; (a)A 7’0 7'7: (b)A€ 7‘0

c’f c’f of \ of

as as C"s Us
7'7: (c)A€M 7'o 7‘2- (d)c 7‘o ”"2: (e)eM 7'0 7‘2- (OM To

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

Figure 3.13: Correlation of effective stress distributions with a type 11 structure map.
The parameters used are: E = 50 GPa, a = 0.05, 0 f = 403 = 300 MPa and ro/rz- = 5.
The effective stress distributions corresponding to various loads in the structure map
are shown as subplots. Points 1—10 represent loads on the four bounding loops,
marking the transition from one phase distribution type to another. Typical phase
distributions are shown for loads at points a—f. Values of (Pi, p0) used for these points
are listed in Table 3.1

65

 

6 .

 

Loading Structure Map
- - - Unloading Structure Map

 

 

 

 

10

o
Pi/cju

Figure 3.14: Structure maps for loading and unloading. Here E = 50 GPa, a = 0.05
and 0RT = 00 = 75MPa, 0.5,RT = 200,05T = 300,0FT = 400, and r0 = 2ri.
The solid curves corraspond to the loading structure map, and the dashed curves
correspond to the unloading structure map. At point 1, p,- = p0 = 0; at point 2,
p,- = 200 and p0 = 400.

66

Chapter 4

Plane Strain Analysis of Shape
Memory Annular Plates Globally

Loaded from A Phase

In this chapter, we will analyze a shape memory annular plate under plane strain
setting following a similar fashion used in Chapter 3. The out-of-plane stress can be
found by the plane strain condition, nevertheless, it can not be expressed explicitly
with the in-plane stresses, neither can effective stress. However, the governing equa-
tions can also be expressed in terms of the effective stress and an auxiliary variable.
Different from the plane stress setting, the Poisson’s ratio affects the stress distri-
bution in this plane strain setting, while the nonincreasing property of the effective

stress along radial direction remains.

67

4. 1 Problem Statement

The problem under consideration is similar to the one in Chapter 3, except that
here we use a plane strain setting instead of a plane stress setting. Therefore, the
out-of-plane stress is not necessary zero, but the out-of-plane strain is zero. Only ax-
isymmetric radial displacement W6") is non-trivial and both azimuthal displacement
and out-of-plane displacement vanish, i. e., 119 = uz = 0. The solution of the plate
involves an, 099, and 032, and they are functions of r only.

The boundary conditions, the equilibrium equation, the strain definition, and the
compatibility condition, are the same as those in Chapter 3. They are included here

for convenience. The boundary conditions are

 

Urr = —po at 7' = To,
(4.1)
Orr = —p,- at r = ft-
The equation of equilibrium is
do 0 — 0
rr + rr 00 = 0. (4.2)
dr r
The infinitesimal strains are
d
err = % and 599 = %. (4.3)

68

The compatibility condition is

 

(1 re
Err = ((136). (4.4)
The plane strain condition requires
Ezz = 0, (4.5)

and this equation is used to find the out-of-plane stress.

From equation (2.12) to (2.14), the stress-strain relation for this plane strain setting

can be written as

 

 

 

f
l a.” : Orr W130 + 022) + 593-,
099 - V(0'7~7~ + 032) t
f 599 = E + Ear, (4'6)
L 522 = Uzz — “(gr + 006) + 5?-

 

Note the transformation strain now also depends on Uzz- The stress tensor, the stress

deviator, and the effective stress for the axisymmetric setting are

69

T = Orr-ref ® 97' + 06690 ® ea + 02292 ® e2)

[(20rr — 0'06 -‘ Uzz)er ® er ‘1' (2060 '— O'zz — Urr)€0 ® 86

OOIH

+(20zz " Urr — Ugolez 8’ 92:],

 

 

\/2
U = 'fflarr "" 066)2 + (0'60 "" Uzz)2 + (022 "" Urr)2.

The transformation strain is obtained by specializing (2.14) to the present case of

plane strain using (4.7) as

tr _ 20m: — 000 - 022

 

 

ET — 20 0&9
59 = 2099 - 3:7, - Orr (,5, (4.8)

 

I 2 — —
( sf," = a” g: 0993115.

Again we restrict the study to the forward transformation starting from pure austenite
and constant temperature at T > A f, and the evolution of martensitic fraction 5 will
also follow equation (3.9).

Use equation (4.6), (4.8), and plane strain condition ezz = 0, 022 can be found as

_ (20'1/ + EOE) (Orr + 0'99)
022 _ 2(0 + Ea5) '

 

(4.9)

Note the right-hand-side of (4.9) contains 0 and 5, which depends on 022, hence this

nonlinear equation does not give an explicitly. As a result, it seems that the effective

70

stress 0 can not be expressed in terms of 0” and 099 directly. Thus it is not obvious
whether one can obtain an explicit expression of d099/dr in terms of arr and 099
as was done in (3.17) for the previous case of plane stress. The next section will use

change of variables to avoid direct evaluation of 0,22 and 0 with given Orr and 099.

4.2 Reformulation of Governing Equations

Introduce the following change of variables:

5rr = 0rr — 0‘22,

(4.10)

' 599 = 099 — 0’22.

Effective stress defined in (4.7) and transformation strain defined in (4.8) can be

rewritten as

 

0 = 5,2... — 07-1099 + 630, (4.11)
and
tr __ 25m" — 569
ET 20 06’
2- _ _
59 = 009 Orr at. (4.12)
20
t 5rr ‘1‘ 5
527‘ = ———2—U-ig-a5.

From equation (4.6) and (4.12), the total out-of-plane strain 822 is

1— 21/ V a __ -
Ezz = E Uzz ‘ (E + 7:)(0rr + 006)- (4'13)

 

71

 

The plane strain assumption ezz = 0 gives

E V a5 _ _
1_ 21/(E + 5X0” + 099), (4.14)

 

022 =

here we assume the range of the Poisson’s ratio is —1 < V < 1 / 2. Strain an and ezz

can be expressed in terms of (Inn, 099 using (4.6), (4.12) and (4.14) as

 

1+V 3a5 _
Err =( ‘l’ )UTT'1

 

E 20
(4.15)
1+V 3a5 _
n..—( E + saline-

The equation of equilibrium (4.2) and compatibility equation (4.4) can be rewritten

as

 

 

 

 

dam. d 1 Ea€ _ _ 5w - 509 _
dr drl1—2V(V+ 20 Harri-099)] + r _0’
(4.16)
1+V 3a5 _ _ d 1+V 3a5_ )_
( E + ’27)(”" 099’ ’3 l( E + 20 >099] ‘ 0‘
Introduce further the following change of variables:
'( 2 7r
l 0 = ——sin + — 0,
7'7. \/'3' (I6 6)
(0 > O) (4.17)

_ 2 , 7r
099 = W Sln(fi — 5W

where 6 is a function of r and for convenience is taken as obeying 0 g 6 < 2n.
This transformation is valid when 0 74 0. When 0 = 0, both 0” and 099 are

necessarily zero, which implies Orr = 099 = ozz. It follows from (4.6) and (4.8) that

72

522 = (1 — 2V)07~r/E = 0, hence 07",- = 099 = 022 = 0. Substitute (4.17) into (4.14)

to obtain

2E V
azz— - T—2_—V(E+2 )0 smfl (4.18)

In terms of 0 = 0(r) and f3 = 5(7), equation (4.16) becomes

f I

 

 

 

 

 

 

 

[(V+E;g)sinfl+1\/_V sin 11([3+%)]E1d-:+
Eag 2” 1r d5 (1—2V)0cosfi
‘7 (V+——)cosfi+——cos(f3+—) ——_—_ ,
4 [ 20 \/3 6 ldr fir (4.19)
3Eg ’ d
(1+u+ 209 )sm m(g_%)§+
3an ii 45 ocosri 3Eag
\H(1+V+ 20900) 8(5—6)$= r (1+V+ 20 ),

where g’ = dg(0)/d0.
Solving the system (4.19) for d0/dr and dfi/dr gives

f

d0 _ _cosz 6K1 — V2)02 + (2 — V)Ea0g + 3E2a2g2/4]
dr rA ’

dfi __ 0cosfi 2 I I -
d; — 4rA {[4(1 V )0 + 3Eag(Eag + 1) + (5 4V)Ea0g )] sm 6+

 

 

\/3(1 — 2V)Ea(g — 0g’) cos 6)},

 

 

(4.20)
where 0 > 0 is required, and A is given by
1 — 2 3E2 2 1 + E
A= V 0+ a gg’+(—L)£(g+09')+
2 8 4
(4.21)

(1 ‘ 24")E“{ [1 + sin(2fi + %)]g + [1 -— 3111(25 + 9109'}-

 

73

The dependence of governing equations (4.20) and (4.21) on parameters E and a is
also solely through the product Ea.

By definition a > O and E, a, g, and g’ are non-negative, also the range of the
Poisson’s ratio is —1 < 1/ < 1/2, and it follows from equation (4.21) that A > 0.

From the first equation of (4.20) it is easily seen that

_ < 0. (4.22)

Thus the maximum effective stress amax occurs at 7‘ = r,- and the minimum effective
stress 0min occurs at r = r0. Accordingly g (1") is also a non-increasing function of
r, i. e., dé/d'r g 0. Now we have similar phase distributions as those in Chapter 3:
six possible phase distributions A, M, (‘3, .46, (3M, ABM will be correlated to load
values p,- and 130. Once again the circle 7‘ = TA in the plate is the boundary between
an A ring and a (‘3 ring, and the circle 7‘ = ”I‘M is the boundary between a C ring and
an M ring.

When the plate is subject to equi-biaxial tractions Pi = p0 = p, the tractions required
for initiating and completing martensitic transformation are very different between
the plane stress situation and the plane strain situation. For plane stress, we have
or,» = 099 = —p and a = |p|. For plane strain, we have arr = 099 = —p, and from

equation (4.9) the out-of-plane stress is

_ —(201/ + Ea§)p
Uzz — 0' + Eaé . (4.23)

 

74

With the third equation of (4.7), the effective stress is given by

(1 - 21/)0

4.24
a + Eafi Ipl’ ( )

a: arr -Uzz| =

and the applied traction p can be retrieved as

0+Ea€
= :t———. .
p 1_ 21/ (4 25)

Therefore, for equi-biaxial traction the pressure 193 to initiate the martensitic fraction

in the plate is

Us

__ 4.26
1 — 21/ ( )

P3 = It
and the pressure pf to complete the martensitic fraction for the whole plate is

0f+Ea 0f+Ea

_ i——— ._ —— . 4.27
pf 1— 2U 03 p8 ( )

Note both 193 and pf are independent of the plate geometry, and the ratio between
them is independent of the Poisson’s ratio. The Poisson’s ratio affects both 193 and
p f, and when the ratio is close to 0.5 it can take very large pressure to initiate and
complete the martensitic transformation. Most of the shape memory materials have
an Ea much larger than of, it follows that p f will be much larger than of. Therefore,
for equi-biaxial loading in a plane strain setting, the pressure required to convert the

whole plate into M phase is much larger than that in a plane stress setting.

75

4.3 A Fully Austenitic Plate (.A Phase Distribu-
tion)

This section considers plates that are completely in the austenite phase and so cor-
respond to a phase distribution of type A. This requires a g as everywhere in the

plate, hence g = g' = O, A = (1 — V2)0'/ 2 and equation (4.20) reduces to

do 2acos2fi d6 2sin 3 cos B
-———— and —=————.

d—r- — 1* dr 'r (4.28)

Note this equation is identical to equation (3.26). The solution to this equation hence

is the same as (3.27), and 57-7 and 599 are given by

for cosfi 31$ 0,

and 6 = c — —Cl
66 1
«$02?2 (4.29)

_ _ C1
Orr—Cl+-\/—§C;T—2’

érr = 699 = :tC3 for cosfl = 0,

where c1, c2 and C3 are constants to be determined by the traction boundary condi-
tions.

For A phase, 6 = 0 and equation (4.14) gives

V

 

 

 

03; = 1_ 21/(51'1' + 660), (4.30)
hence
2
azz = V61 for cosfi aé 0,
1 — 21/
(4.31)
2
022 = W3 for cosfi = O.
1 — 21/

76

 

Substitute (4.29) and (4.31) into (4.10), radial stress arr and hoop stress 099 are as

 

 

 

 

follows
Cl Cl Cl C1
0 = — + and 0 = — for cos 0,
"T 1.. 2,, «302.2 99 1- 2U @212 B 7‘
arr = 099 = :1:1 :321/ for cosfi = 0.
(4.32)

If p,- # 190, then the constants c1 and oz can be obtained from the first equation of

(4.32) as
(1" 2V)(Pi7"2 — 1907(2))

C —- 1
1 .342

 

for 191 ,1 P0 (433)
= (1" 2V)(pi7'i2 _ p070)

73(190- Fargo?

 

If instead p,- = p0 = p, then the second equation of (4.32) gives C3 = (1 — 21/)|p|

whereupon

0rr(r) = 090(7) = -P. 0(r) = (1 - ZVllpl for Pi = 190. (4-34)

Equations (4.32), (4.33) and (4.34) together retrieve the classical linear elastic solu-

tion, with stresses given by

 

 

 

 

 

,
a _ pin-2 - par?) (po- war???»
N 7% — r12 (7‘3 — r?)r2 ,
< 0 _ 1327.12 _ P073 _ (pO- p2)r2rg (4 35)
00 7'3 -— r22 (r3 — 7;?)1‘2 ,
21/(pirz-2 — porg)
022 = 2 2 -
1 T0 - Ti

 

77

The effective stress is given by

 

A \/(1— 21/120322. -r,?p.>2+3r;1r3<p.- -po)2/r4
a = 0(r) := 2 2 . (4.36)

 

The maximum and minimum effective stresses are given by amax = 6(7)) and 0min =
6(7‘0) respectively.
The validity of this solution requires omax S as. The case amax = as is associated

with parameter values obeying

[(1 — 21/)2rf‘ + 37gb)? — 2rg[(1 — 21/)27‘7;2 + 3rg]pipo+
(4.37)

4(1 — V + V2)7“01pg — (1"(2) — 1",?)203.

Equation (4.37) defines an ellipse in the (pi,po) plane named FAL as described in
Chapter 3. The FAL is symmetric with respect to 180° rotation about the origin, a
consequence of the tension / compression symmetry in the modeling description.

An alternative way of obtaining the above FAL is to solve (4.29) subject to a require-

ment 0(r,) = as. If B(r,-) = 60 7g 7r/2,37r/2, then according to (4.29) cl and C2 are

78

given by cl = as sin 60 and c2 = tan 60/732, hence the stresses are given by

2

cosfior-
6 =0 sin +———7’,
rr 3( 30 \/37'2 l

_ . cosfiori2
000 = 03(Slnfi0 — 7372—),

 

4
r.
a = as\/sin2 fig + 0032 60—2,
< T (4.38)

 

a -— sin 3
a ,
2:2 1 2 S 0

sinfio cos 50732
0rr=03(1—_—2u+ 73372 —)’

a _0(sinflo cosfiorzz)
( 66— 81—21! «3.72 .

 

Note this solution is also valid for 60 = 7r/ 2,37r/ 2. Applied pressure p,- and p0 are

then given by

 

_(_ sin 60 cos 50 sin 60 COS flag-2

1 _ 21/ J3 )as and p0 = -( 1 _ 2V fira 2..)0 (4.39)

 

:91:

Equation (4.39) defines an ellipse in parametric form within the (19,3190) plane where

,60 is the varying parameter. Solving (4.39) gives sin ,30 and cos 50 as

(1 — 2140947"? — P073)

 

 

Sinfio= 1
(7221 — rizlas (4 40)
./3( — -)r2 '
cosfl — po 1?, 0
0“ 2_ 2 ’
(7'0 Tilas

whereupon the identity sin2 fig + cos2 60 = 1 retrieves (4.37).

79

 

The following will use fi as a varying parameter in obtaining loops associated with
the remaining transitional cases TA = r0, ”PM = 7‘2', TM = To to produce a complete
structure map. The structure maps and stress distributions are similar to those for

the plane stress situation.

4.4 A Fully Martensitic Plate (M Phase Distribu-
tion)

This section considers plates that are completely in the martensite phase and so
correspond to a phase distribution of type M . Hence 9 = 1 and g’ = O, and equation

(4.20) reduces to

 

do _ _c052 fi[(1 — 112)02 + (2 — V)Eaa + 3E2a2/4]

 

 

 

E: _ TA
(4.41)
dfi _ acosfi 2 ,
d7 — 47A [4(1— 11 )0 + 3Ea] s1nfi + \/§(1— 2V)Eacos fi)},
where a > 0 is required, and A is given by
_ 2 _ _
A = 1 2V 0 + (2 :)Ea + (1 241/)Ea sin(2fi + %)' (4-42)

The equality case for the validity condition 0(7'0) Z of defines the FML in the (191,190)
plane. The FML at fixed of, Ea, V and ro/ri can be obtained from integrating
equation (4.41) with boundary condition (0(ro), fi(ro)) = (of, fig) upon varying fig

from O to 27r. Like the FAL, the FML is also symmetric with respect to 1800 rotation

80

about the origin. The FML corresponds to the transitional condition TM = 7'0, and
the region exterior to the FML corresponds to tractions (1923290) that generate an
M phase distribution. The FML will surround the FAL in the (1923190) plane, and
the FML is generally not an ellipse. The region between the FAL and the FML
corresponds to loads (1723170) that generate one of the phase distributions containing
a mixture zone: (3, AG, CM or ABM.

Although the FML is usually determined numerically, it will pass through (12,3130) =
(p f, p f), where p f is defined in (4.27), which provides the uniform equi-biaxial stress
solutions on the F ML. More generally, the line p,- = p0 provides the uniform equi-
biaxial stress solutions in the (pi, p0) plane. On this line we have (Try-(7‘) = 099(r) =
—p, and a uniform effective stress is given implicitly by a nonlinear equation (4.25).
The phase fraction field {(7') is also independent of r and can be determined after a
is found. Such loading results in a phase distribution of either type A, C or M.
Figure 4.1 displays both the FAL and the FML for the parameters E = 50 GPa,
a = 0.05, V = 0.3, as = 75 MPa, of = 300 MPa, and ro/ri = 2. Since the FML is
much larger than the FAL, Figure 4.1(b) is the enhanced view of the zone enclosed in
the dashed rectangle in Figure 4.1(a). Note also that the FML is highly elongated in
the p,- = p0 direction, which suggests that in order to produce an M phase distribution
efficiently, we need to avoid a biaxial loading situation. Also 19,; and p0 appear to be
equally effective in producing an M plate, which is quite different from the plane
stress situation where p0 is more efficient than Pi in producing an M plate.

The following will study the effects of material parameters on the FML. An F ML is

completely determined by the parameters of, Ea, V and ro/ri. Figure 4.2 shows the

81

 

 

 

 

 

 

 

 

 

100'.'I'I'I'I'I'l'lll'4 20 'fi ' I 1
80 - 4 _/
4o '— l '0 ' ‘
20 l - J
po/O's 0L 4 po/O'sfl- / -
'20 .' 1 i
‘40 T '- 40 -
-60 - _ F
_80} _ . /
_100 A l A l A l A l A l A l A l A I Ll A 4 -20 l I l A
-100-80-60 .40-20 0 20 4o 60 80100 -2o -10 o 10 20
P7; /°'s pi /0’3
(8) (b)

Figure 4.1: The FML (Fully Martensitic Loop) and FAL (Fully Austenitic Loop) for
E = 50 GPa, a = 0.05, V = 0.3, of = 403 = 300 MPa and ro/rz- = 2. Tractions inside
the smaller ellipse (FAL) produce a fully austenitic plate, while tractions outside the
bigger loop (FML) produce a fully martensitic plate. Part (b) is the enhanced view
of the zone enclosed in the dashed rectangle in (a).

effect of varying of on the F ML at fixed Ea, V, and ro/ri. In particular, this figure
displays a family of four F MLs corresponding to O'f/O's = 1, 4, 8, 12, with the other
parameters set at E = 50 GPa, a = 0.05, V = 0.3, as = 75 MPa and uh; = 2. As
one would anticipate, the FMLS are nested with respect to increasing of. Note for the
plane stress situation, the FML corresponding to O'f/O's = 1 intersects the FAL at the
equi-biaxial stress solutions (pi, p0) = :l:(af, of). But this does not hold for the plane
strain situation, as it is easily seen from equation (4.27) that p f = (1 + Ea/03)p3 if
0f = as.

Figure 4.3 shows the effect of varying Ea on the F ML at fixed V, of and ro/ri. In
particular, this figure displays a family of five FMLS corresponding to Ea / 500 MPa
= 1, 2, 5, 8, 10, with the other parameters set at V = 0.3, 0f 2 403 = 300 MPa,
and ro/rz- = 2. An FML for a larger Ea encloses an FML for a smaller Ea. When
Ea = 0, equation (4.41) reduces to (4.28) and FML is an ellipse but different from

82

120

 

 

 

 

_120.1.1L111411
-120 -80 -40 0 4O 80 120

pi/O's

 

Figure 4.2: FMLS for various of values. Here E = 50 GPa, a = 0.05, V = 0.3,
as = 75 MPa and ro/rz- = 2. All the F MLs are marked with corresponding value of

Uf/Us.

FAL is of > as. We do not plot this FML here because it is very small compared to
the other F MLs.

Figure 4.4 shows the effect of varying V on the F ML at fixed Ea, of and ro/ri. In
particular, this figure displays a family of five FMLS corresponding to V = O, 0.05, 0.1,
0.2, 0.3 with the other parameters set at E = 50 GPa, a = 0.05, of = 403 = 300 MPa,
and ro/rz- = 2. An FML for a larger V encloses an FML for a smaller V.

Figure 4.5 shows the effect of varying ro/ri on the FML at fixed Ea, V, (If and Ea.
In particular, this figure displays a family of five F MLS corresponding to ro/rz- = 1.2,
1.6, 2.0, 2.5, 3.0, with the other parameters set at E = 50 GPa, a = 0.05, V = 0.3,
and of = 403 = 300 MPa. All the five FMLS pass through the two common uniform
equi-biaxial load points (pi,po) = :t(pf,pf) = (93.303,93.303) (recall from (4.27)

that p f is independent of 7‘). For small ratio of ro/ri, for example, ro/rz- < 2.5, the

83

 

 

ImLAJ

1 1‘4

4'1

 

1
.1

 

 

l A l A l A l

00 ‘
-200 -150 -100 -50 0
191/03

A l A l A l

50 100 150 200

 

Figure 4.3: FMLS for various Ea values. Here V = 0.3, of = 403 = 300 MPa and
ro/rz- = 2. All the FMLS are marked with the value of Ea/ (500 MPa).

FMLS elongate significantly along the equi—biaxial direction. When increasing ro/ri,

the FMLS expand dramatically in p,- direction.

4.5 Stress Distributions and Structure Maps

We now add two other loops to complete the structure map: ADL, which is associated
with 0(ro) = as, and MEL, which is associated with 0(r,) = (If. With the exception
of locations on the equi-biaxial line, both the ADL and the MEL are determined
numerically. This is most easily accomplished using the same procedure that generates
the FML upon replacing the FML condition 007;) = of with either the ADL condition
0(ro) = 03 or the MEL condition 0(r,) = of. The ADL and the MEL are generally
not elliptical, and they are symmetric with respect to 1800 rotation about the origin,

hence the symmetry of structure map follows.

84

[00

 

50

pO/OS 0

-50

 

 

 

 

pi /Gs

Figure 4.4: FMLS for various V values. Here E = 50 GPa, a = 0.05, of = 403 =
300 MPa and 1‘0 / r,- = 2. All the FMLS are marked with the value of V.

Two examples of fully articulated structure maps are presented in Figure 4.6 and
Figure 4.7. Since F ML is much larger than FAL, Figure 4.6, 4.7 and 4.8 have either
two or three parts. In each figure, part (a) is the full view of the structure map; part
(b) is the enhanced view of the zone enclosed in the dashed rectangle in part(a); part
(c),if it exists, is the enhanced view of the zone enclosed in the dashed rectangle in
(b)-

Both Figure 4.6 and Figure 4.7 correspond to material parameters E = 50 GPa,
a = 0.05, V = 0.3, as = 75 MPa, of = 300 MPa. Their difference is due to differing
values of ro/ri, namely ro/rz- = 2 for Figure 4.6 and ro/ri = 8 for Figure 4.7.
Each structure map is partitioned into various regions by the four loops FAL, FML,
ADL and MEL. Each region corresponds to a distinct phase distribution type as also
displayed in Figure 4.6 and Figure 4.7.

The qualitative difference between these two structure maps is that in Figure 4.6

the ADL is strictly inside the MEL while in Figure 4.7 it is not. We shall refer

85

 

100 ' I ' l ' I ' I ' I

 

 

 

 

50 ~ .
130/03 0 l" _,
J

-50 _ _

l
-100 A 1 A 41 A LA LA 1 AL LL, A
-150 -100 -50 0 50 100 150
1972/03

Figure 4.5: FMLS for various ro/ri values. Here E = 50 GPa, a = 0.05, V = 0.3 and
of = 403 = 300 MPa. All the FMLS are marked with ro/rz- values.

to the former as a type I structure map and to the latter as a type II structure
map. A type II structure map has two symmetric regions corresponding to an ABM
phase distribution whereas a type I structure map does not involve this type of phase
distribution. To within unit scaling, the structure map is determined by the material
parameter combinations Ea, V, of, 03 and geometry parameter ratio ro/ri. For fixed
material parameters Ea, V, of, 03 one finds that the structure map is of type I for
small ro/rz- but is of type II for large ro/ri. The condition To sufficiently greater than
r,- makes possible a situation supporting the condition 7°,- < TM < TA < To that is
associated with the phase distribution ABM.

For given Ea and of, as there is a special value of ro/ro such that the structure
map is of type I below the special value but of type 11 above the special value. For
the material parameters associated with Figures 4.6 and 4.7 (Ea = 2.5 GPa, V = 0.3,

of = 403 = 300 MPa), this special transitional value is found to be ro/ri = 5.52. The

86

 

 

 

 

 

 

 

 

 

 

100 IO
50 h 5 y.
170/630: po/an-
-50 -5 ..
_loo . A A A A I A A A A l A A A A I A A A A ll -10 A
~100 -50 0 50 100 -10 -5 0 5 10
pi/Os pi/OS
(a) (b)

Figure 4.6: Structure map (type I) for ro/rz- = 2. Here E = 50 GPa, a = 0.05,
V = 0.3 and of = 403 = 300 MPa. Part (b) is the enhanced view of the zone enclosed

in the dashed rectangle in (a). Note that there is no ABM region. Traction pair
(19,3170) = (8, 1)03, marked as a dot, produces a BM phase distribution. The resulting
stresses are shown in Figure 4.9.

structure map for Ea = 2.5GPa, V = 0.3, of = 403 = 300 MPa, and ro/ri = 5.52 is
plotted in Figure 4.8 showing how the ADL osculates the MEL.

The stress distribution for the phase distributions under discussion can be determined
numerically on the basis of (4.20) using the iterated shooting method. For example,
Figure 4.9 displays the radial variation of a, arr, 099 for the BM phase distribution
associated with the load (pi,po) = (8,1)03 represented by the dot in Figure 4.6
(Ea = 2.5 GPa, V = 0.3, of = 403 = 300 MPa, ro/rz- = 2). A mixture of austenite
and martensite appears on the outer plate boundary 7" = To. This mixture becomes
progressively richer in martensite as 7' decreases to r = TM = 1.3613- within which
the plate is in the martensite phase. This solution could, in principle, be obtained
by many quasi-static load paths in the (1923190) plane so long as the effective stress

at each point in the plate remains non-decreasing. In particular, the proportional

87

 

 

 

 

 

 

 

 

100 -

 

 

6°
\

«9
O
I

 

 

  
  

M

I L I A I A I A I

-100 -50 0 50 100 ~20 -10 0 10 20

 

 

 

 

 

 

Pi /63 pi /0.9
(b) (C)

Figure 4.7: Structure map (type II) for ro/ri = 8. Here E = 50 GPa, a = 0.05 ,
V = 0.3 and of = 403 = 300 MPa. Part (b) is the enhanced view of the zone enclosed

in the dashed rectangle in (a); part (c) is the enhanced view of the zone enclosed in
the dashed rectangle in (b). All six types of phase distributions are present. Traction
pair (pi,po) = (8,1)03, marked as a dot, produces a ABM phase distribution. The
resulting stresses are shown in Figure 4.10.

loading path (Pi: p0) = (8, 1)03k where k is a load path parameter is associated with
an initial phase distribution of type A. There will then be three values of k < 1, say
k1 < k2 < k3, associated respectively with transitions through the sequence of phase
distributiOns: A, AB, B, BM. Here I: = 1:1 is associated with the emergence of the
AB front 7' A at the inner plate boundary. This front proceeds through the plate as

k increases until it merges with the outer boundary when k = kg. The BM front TM

88

then emerges at the inner plate boundary when k = k3 and continued loading drives
this front to r = 1.361'2- when k = 1. Further loading will then drive this front to the
outer boundary, resulting in a phase distribution of type M.

Figure 4.10 shows the stress distributions at the same load (Pi, p0) = (8,1)03 as in
Figure 4.9 but for the structure map of Figure 4.7. In this case the larger value of
r0 (r0 = 873-) gives rise to a phase distribution of type ABM with r A = 6967‘,- and
TM = 1.2513. In this case the proportional loading path (19,3190) = (8,1)03k passes
through phase distribution states of type A, AB, ABM with the transitions taking
place at two special k values that are analogous to k1 and k3 of the previous example
associated with Figure 4.9.

A final set of stress distributions is shown in Figure 4.11. These correspond to the
osculation point of ADL with MEL at (pi,p0) = (4.43,0.68)03. Along the propor-
tional loading line from the origin through this point, this solution is associated with
the disappearance of the pure austenite phase at the outer plate edge (7" A = To) that

is simultaneous with the appearance of the pure martensite phase at the inner plate

edge (TM = Ti)-

89

 

150 ' V v ' 1 Y T V fir I ' V I U 1* fi'

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

 

l A
100 - A
,
50 h M
170/03 0 r
-50 -
-100 _.
-150 h A A A A 1 A A A A L A A A A 1 A A A L l
-500 -250 0 250 500
pi/O's
(a)
100 10 f T A I T I
50
po/O's 0
-50
_100 A l A 1 A l A
400 -50 0 50 100
pi/Gs
(b)

 

Figure 4.8: Structure map for the special value ro/rz- = 5.52. Here E = 50 GPa,
a = 0.05, V = 0.3 and of = 403 = 300 MPa. Lower ro/rz- values produce type I
structure maps, while higher ro/rz- values produces type II. Part (b) is the enhanced
view of the zone enclosed in the dashed rectangle in (a); part (c) is the enhanced
view of the zone enclosed in the dashed rectangle in (b). There is no ABM region,
and ADL contacts MEL at two points. One of the two points (shown as a dot)
represents (pi, p0) = (4.43, O.68)03, producing a B phase distribution with 0(r,) = of
and 0(r0) = 03. The resulting stresses are shown in Figure 4.11.

90

 

 

 

 

 

1.0 1.2
rr

- 6 L+::MF=’m"3o
/

-3 _ ’ / 0“,,

Figure 4.9: Stresses (normalized by as) for a BM phase distribution. Here E

50 GPa, a = 0.05, V = 0.3, Of = 403 = 300 MPa, ro/rz- = 2 and (pi,po) = (8,1)03.
An elastic solution (a = 0) is shown for comparison (dashed lines).

 

 

 

 

 

 

' 20 or = 0.05 a: 0
16
12
8 ‘ \ TM = 1.2572,;
4 _l 590
TA = 6.9613;
0 ‘ :‘- r/rz
2,,~ —:====6 8
/
.4 ..
67'7“
22
-8 -+

 

Figure 4.10: Stresses (normalized by as) for a ABM phase distribution. Here E =
50 GPa, a = 0.05, V = 0.3, of = 403 = 300 MPa, ro/rz- = 8 and (1223190) = (8,1)03.
An elastic solution (a = 0) is shown for comparison (dashed lines).

91

 

 

 

 

 

 

 

 

 

 

Figure 4.11: Stresses (normalized by as) for a special B phase distribution with
TM = r,- and TA = To, for E = 50 GPa, a = 0.05, V = 0.3, of = 403 = 300 MPa,
ro/ri = 5.52 and (pbpo) = (4.43,O.68)03. An elastic solution (a = 0) is shown for
comparison (dashed lines).

92

Chapter 5

Plane Stress Analysis of Shape
Memory Annular Plates Subject to

Arbitrary Traction History

In this chapter, we will study the responses of shape memory plates subject to arbi-
trary loading history in the plane stress setting. Partial martensitic transformation
associated with the transformation subloops may take place. The non-increasing
property of the effective stress along radial direction does not necessary hold, and
the structure maps described in Chapter 3 can not be used for determining the phase
distribution for given traction loads. Numerical approaches are adopted as an an-
alytical study seems unlikely. The given traction history is discretized into a series
of time steps and the numerical simulation is performed step by step. First, the

shooting method is used to solve the responses of an annular plate under a certain

93

class of traction histories, however, this scheme seems to have difficulties in handling
the situation when both loading and unloading occur in the plate within a time step.
Alternatively, an iterative finite difference procedure is shown to be efficient in solving
the proposed boundary value problems. The numerical results showing the history
dependence in stress responses, the non-monotonicity in effective stress along the ra-
dial direction, the repeating or stabilizing responses of the plates under cyclic loading,

and the movement of B /A and M/ B interfaces are given in detail.

5. 1 Problem Statement

For a general traction history, the martensitic transformation can follow any path
as described in the hysteresis model of Chapter 2. Furthermore, the martensitic
transformation can follow different subloops for different locations in the plate. An
analytical investigation on the properties of the effective stress, such as those carried
out in Chapter 3 and Chapter 4, seems unlikely, therefore numerical studies will be
performed.

To aid the numerical simulation, the applied traction history is represented by discrete
time steps. Here we denote the total number of time steps as K, and t k is the time for
the kth time step where 0 _<_ k S K. Note tmm _<_ tk g tmax, tg = tmz-n, t K = tmax,

and t,- < tj if 0 S i <j g K. From equation (3.1), the tractions at t = tk is

94

and ASA-(to) = 15000) = 0.

Although the martensitic evolution function used here is different from the one used

in Chapter 3, the boundary conditions, the equilibrium equation, the strain defini-

tion, the compatibility condition, the stress-strain relation, and the effective stress

definition, are the same for both cases and are included here for quick reference. The

boundary conditions are

07‘7‘ = -p0 at 7' = T0)

07'1' = —pi at 7' = Ti.

The equations of equilibrium is

darr + 07‘1" — 009 =

 

 

0.
dr 1"
The strains are defined by
5 _dur and 5 _ur
rr - d?‘ 60 - 7. '
The compatibility condition is
E _ (10660)
rr — -
(17'

95

(5.4)

(5.5)

The total strains are given by

EM = ___Urr 21/066 + 55",
— V

—V a +0
522 ___ ( 7‘71“; 66) +5272

 

 

where the transformation strains are

2 _.
a": Orr ”660,5,

 

r 20

2 _
531‘ = 009 07'7" a6, (5.7)

20

o + a
Etzr = ...—Wag
The effective stress a is given by

0' = 072»,- — arragg + 036. (5.8)

Since 5 and a are not uniformly distributed in general, it’s possible that partial
transformation will introduce a situation that for different 1" the evolution of trans-
formation follows different trajectories in g — a plane. Assume that the time steps
given are small, and the effective stress at given 7" within a time step can be treated as
either nonincreasing or nondecreasing. Taking those time steps as the time sequence
described in Section 2.1, equation (2.8) and (2.10) can be used to describe the trans-

formation behavior for any 1‘. The equilibrium equation is solved step by step for the

96

given traction history, with the stress and martensitic fraction from the previous step
used for determining the martensitic evolution for the present load step.

Only isothermal conditions will be considered, and we will drop the temperature
dependence of transformation starting and finishing stresses, therefore we omit the
temperature dependence of QFT and QRT on temperature T. For convenience, we
combine (2.8) and (2.10) as g(a,gs) to describe the transformation behavior. The

martensitic volume fraction g is then given by

E = 10.798) (59)

The term gs is load history dependent, but remains as a constant for any time step
since we assume that the load type (loading/ unloading) does not change within any
time step. Equation (2.9) and (2.11) are used for computing gs according to the
loading type at given 7'. Note gs at given 1" is uniquely determined given a and E of
the previous step and the loading type of this step, while the actual value of effective
stress does not affect gs as long as the loading type is determined. In general, gs is
a function of T, which needs attention when taking derivative of g(a, gs) about 7*. If
a plate is globally loaded from the austenite phase, then 9,, E 0; if a plate is globally
unloaded from the martensite phase, then gs E 1.

The compatibility equation in terms of stress can be expressed as

d 20 — a d
g(arr + 090) - WEEM’QS) = 0, (5.10)

97

where E3(0, gs) is the secant modulus given by

1
E3995) = 1/E+aa(a.gs>/o‘ ‘5'“)

 

By chain rule, the derivative dEs(0, gs)/dr is

dEs 6E3 FE,- (9E3 993

ar=aadr+a73dr° (“2)

Substitute equation (5.12) into (5.10) and remove dorr/dr using (5.3), dagg/dr can

be solved and finally we have the following standard form of a two-variable ODE

 

 

 

 

system:
dUrr = ___Ur'r " 099
dr 1' ’
(5.13)
C1066 = B (Grudge) Orr - 096 + C(Umaoo) dgs
dr A(0rr,0'gg) r A(am~,099) dr ’
where
A(a7~7~, 099) = 403 + Ea[3a,gr§ + 06,007” -— 2099)2],
B(0rr, 099) = 403 + Ea[402§ + (g — 06,0)(202 — Barrage”, (5-14)

89(0198)

C(0rr,060) = 2EC¥U2(U7~T — 2006) 89
8

I

and g = g(a, gs), 51,0 = 6§(a,gs)/80. The differential equation (5.13) is valid only
when a, arr, 099, and gs are smooth so that they have derivatives about r. Note

A(a,~,~, 099) > 0 as long as a aé 0, which is obvious since 6 2 0 and Q0 2 0.

98

In our study, the applied normal tractions are taken to be a continuous function of
time and no rate effect is taken into account. Therefore, the corresponding martensitic
transformation in the plate will proceed continuously about time and the martensitic
volume fraction will be continuous about 7' in view of the continuity in the given
martensitic evolution function. In 7“ direction, this will require a continuous a field,
and the continuity of 099 follows from the continuity of arr, which is required so as
to maintain equilibrium.

The history data is known implicitly, which will be used for evaluating 93(7) With
the continuity in 0 and g, the function gs(r) will be continuous and smooth except
for a finite set of r values due to the coexistence of loading and unloading. We use
H to denote the set of points at which one of the derivatives of a, arr, 099, and gs
about 1" does not exists. For a plate globally loaded from austenite phase, or globally
unloaded from martensite phase, IHI = 0 since gs is uniform and a, an», and 099 are
smooth as can be seen from the stress solutions in Chapter 3. In fact, whenever a
plate is fully martensite, we can set gs = 1 and H = 0, and discard all the previous
history dependent data; similarly, whenever a plate is fully austenite, we can set
gs = 0 and H = 0, and discard all the previous history dependent data. Therefore,
equation (5.13) is valid for 7‘, < r < 7‘0 and r 3 11-11.

With the continuity of arr and 099, we can essentially integrate equation (5.13) if
an and 099 are given at one location. If arr and 099 are given, then 0 and load
type can be inferred, and gs can be obtained. The next section will use the shooting
method, which requires the integration of (5.13), to obtain solutions for a specific

class of traction histories.

99

5.2 Numerical Solution by Shooting Method

In Chapter 3, we used the shooting method for solving the stress distributions. Here
again we use this method to solve the plate responses for a traction history producing
global loading solution or global unloading solution, but without the requirement
that the initial phase is austenite or martensite. Recall g3 can be obtained just by
the loading type, then the gs distribution can be computed with the given loading
type for the whole plate without knowing the actual effective stress distribution.
The derivative dgs/dr can be approximated by finite difference scheme, and equation
(5.13) can be integrated. Different from the analysis in Chapter 3 where dgs/dr = 0,
here gs is not necessarily constant therefore its derivative about 7” is not always zero.
A numerical example is given here for a plate with ro/rz- = 2, E = 50 GPa, a = 0.05,
afiT = 00 = 100MPa, 05,71 = 200, UfT = 300, and 0;? = 400. A traction history
under which the plate undergoes global loading for 0 S t S 1 and global unloading
solution for 1 S t S 2 is given by

151'“) = 5tag and 130(t) = 3tag for 0 S t S 1,
(5.15)

15,;(t) = 5(2 — t)ag and 130(t) = 3(2 — t)ag for 1 < t g 2.
Figure 5.1 shows this path on the loading structure map for the annular plate, and a
total of 40 time steps are used, with 20 for loading (0 _<_ t g 1) and the other 20 for
unloading (1 _<_ t S 2). Although more time steps can be used, we use a step number
of 40 for clearer plotting. The global loading and global unloading conditions for this

traction history are also verified for, larger step numbers. There are 21 points on the

100

 

 

 

 

 

-6 -21 ' $2

Figure 5.1: Load steps in (pbpo) plane. Here E = 50 GPa, a = 0.05, of” = 0'0 =

100MPa, GET = 200, afT = 300, of? = 400, and r0 = 213;. The curves on the

plate are from the loading structure map. A total of 40 time steps (At = 0.05) is
used. For 0 < t S 1 (step 1 to step 20), the plate is globally loaded to ABM phase
combination. For 1 < t 5 2 (step 21 to step 40), the annulus is globally unloaded to
A phase.

figure as the loading and unloading paths are overlapping. At t = 1, (pi, p0) = (5,3),
which locates inside the ABM region of the loading structure map. The ensuing
unloading will introduce partial reverse transformation since there is B phase in the
plate.

Figure 5.2 shows the effective stress (normalized by 00) for all the time steps. In
Figure 5.2(a), from bottom to top, the effective stress curve corresponds to initial
step (t = 0) to step 20 (t = 1). For this time period, the plate undergoes global
loading since throughout the plate the effective stress is increasing. In Figure 5.2(b),

from top to bottom, the effective stress curve corresponds to step 20 (t = 1) to step 40

101

 

 

 

 

 

 

 

 

 

 

(a) . . . . (b)

Figure 5.2: Effective stresses (normalized by 00) by the shooting method for the
loading steps correspondin to Figure 5.1. A total of 40 time step is used, and
E = 50 GPa, a = 0.05, of = 00 = 100MPa, 053T = 200, of T = 300, of T = 400,
and r0 = 21",. (a) Effective stresses for 0 S t S 1 (initial step to step 20), the plate
undergoes global loading. (b) Effective stresses for 1 S t S 2 (step 20 to step 40), the
plate undergoes global loading.

(t = 2). For this time period, the plate undergoes global unloading since throughout
the plate the effective stress is decreasing.

Figure 5.3 shows 6— 0 plot for 0 S t S 2 at three locations, 1" = 1, r = 1.5 and r = 2.
At T = 1 (circle points), it first undergoes forward transformation from austenite to
martensite, then it is transformed back to austenite; at r = 1.5 (diamond points),
it first undergoes forward transformation from austenite to a martensite/austenite
mixture, then it is transformed back to austenite; at r = 2 (star points), there is no
martensitic transformation taking place. This plot illustrates a scenario that different
locations can follow different trajectories on 5 — a plane.

Theoretically, the shooting method should also be applicable for arbitrary loading
history. Although significant amount of effort has been devoted to developing various

schemes based on directly integrating equation (5.13), no robust algorithm has been

102

1.0- ’3

    
   

°<><><><><>0<><>o

 

 

0' ‘ 100'. l 200001 360 460

Figure 5.3: Evolution of martensitic fractions for 1' = 1, r = 1.5, and 1' = 2. Effective
stress a is normalized by 0'0. The applied load steps correspond to Figure 5.1. A

total of 40 time step is used, and E = 50 GPa, a = 0.05, of” = 0’0 = 100MPa,

of” = 200, JET = 300, of? = 400, and 1'0 = 21“,.

found. The primary difficulty comes from determining the loading type as the inte-
gration proceeds along radial direction, especially when the effective stresses for two
consecutive time steps are very close. Since the loading type is essential in evaluating
gs, dgs /d1~ and g, a false loading type due to numerical error can invalidate the results
of an integration. Finite difference method, such as relaxation method (Press at al.,
2002) and iterated deferred correction (Cash, 1986, Daele and Cash, 2000) are can-
didates to solve (5.13) by iteration, however the evaluation dgs / (11‘ needs numerical
approximation, which could be problematic when two mesh points have different load-
ing type. Alternatively, an iterative finite difference scheme that avoids the evaluation

of dgs/dr is introduced in the next section.

103

5.3 Numerical Solution by Iterative Finite Differ-

ence Method

Given a transformation strain field, the stress field can be integrated from the govern-
ing equations with the associated integral constants determined from the boundary
conditions. The task is then to find the transformation strain field that produces a
stress field (based on the integration) such that these two fields satisfy the constitu-
tive relation exactly. A nonlinear equation on transformation strain is obtained, and
it is discretized by a finite difference scheme and solved by numerical iterations. This
scheme is robust and the numerical results from this scheme agree with the available

results from the shooting method.

5.3.1 Stress Solutions by Integration

Substituting equation (5.6) into (5.5), we have
- V0 0 — Va
—-——Urr 60 + eff = {r[——09 E W + E?]}'. (5.16)

E

This equation can rearrange to
1 I I t'r ’ tr
E[(1 + V)(am~ — 099) - r009 + wow] = (1‘56 ) — 5,. . (5.17)

With equilibrium (5.3), the V in (5.17) vanishes and the equation is simplified to

104

1 I
g(t)” -— 099 — r089) = (r5?) - 55.7".

From equation (5.3), 099 is

_ I

Substitute (5.19) into (5.18), we have an ODE on arr

r2042. + 3171;, = E [Eff — (TEETH,

and this equation can be rewritten as

Integrate equation (5.21) once, we have

130;". = H(1°) + 61,

where

Hm = E1] 915%) + sword.) + 7.263%.) — am},

and 61 is a constant. Immediately,

0;. = [Hm + 611/23.

105

(5.18)

(5.19)

(5.20)

(5.21)

(5.22)

(5.23)

(5.24)

Further integrating (5.24) gives

a... = /. H(9)/93d9 — eA/(2r21 + co.

Finally, we have or,- as

arr = Cg + 61/7‘2 + [(1"),

where

1(7) = [mm/93d.)

and C1 = —Cl /2. Hoop stress 006 can be obtained via (5.19) as
099 = Cg + [H(T) — C1]/1'2 + [(7‘).
From equation (5.26), the boundary condition (5.2) becomes

Cg + C1/1‘z.2 = —p,',

Cg + C1/r3+ 1(1‘0) = —p0.
Solving (5.29) for constants Cg and C1 gives

_ P1732 — [Po + [(70)ng
— 2 2 ’
[Po - 191+ 100)]?ng

C =
1 .g_.,2

 

 

106

(5.25)

(5.26)

(5.27)

(5.28)

(5.29)

(5.30)

Note for the above deduction of the stress solution, we only require that the given
transformation strain is integrable, which is a fairly lax condition. With the continuity
of martensitic volume fraction and stresses, the transformation field, calculated by
equation (5.7), will be integrable because there is no singularity involved.

If there is no transformation strain, we have H (1') = I (r) = 0, hence

 

 

171'"? ‘19ng
CO = 2 2 1
2 2 (5.31)
C _ (Po-Pilriro
1 .3 -..2 '

Z

The classic elastic solution is immediately retrieved by substituting (5.31) into (5.26)

 

 

 

 

and (5.28) as
pin-2 - 1201-3 (P0 - Film-27%
0"” = 2 2 + 2 2 2 ’
r0 — r2. (1'0 — 1'2. )1' 5 32
.72 _ 7.2 ( _ _)7.27.2 ( . )
099=pzi p00_ P0 131 1' 0.
1‘3 — 1'22 (1% — 1‘i2)1‘2

Given pi, p0 and any estimation for the transformation strains, the radial and hoop
stresses can be found as equation (5.26) and (5.28) respectively. With these stresses,
the martensitic fraction field 5 can be computed by (5.8) and (5.9). Note (5.9) requires
to determine to invoke (2.8) or (2.10) at each 1‘. After the 5 is computed, one may now
retrieve transformation strains from (5.7). If this computed transformation matches
the estimation, then we found a solution for the given traction boundary condition.
Otherwise, this computed transformation strain can serve as an improved estimation

for the transformation strains.

107

During this procedure, integral H (1‘) and 1(1‘), constants Cg and c1, stresses on,
099 and a, martensitic fraction 5, and (new) transformation strains eff and 53" are
uniquely determined by pi, p0 and the previously given transformation strains. We

tr tr tr

denote this procedure as two functionals, Er(5.,. ’56 ,1“) to compute the (new) 8,.

field, and E6 E”, 5”, 1" to compute the new Etr field, that is
7' 9 0

613‘ = E703". 53", 7"),

(5.33)

5th = E6053", 53", r).

In general the transformation strain distributions are not known beforehandoexcept
for a fully austenite plate where the transformation strains are trivial. An issue
remains is how to select between (2.8) and (2.10) in the above algorithm, otherwise
the equation (5.33) may not be uniquely determined. This is to be expect as it is a
direct consequence of the path dependence of the hysteresis algorithm (2.8) and (2.10).
A convenient way to resolve this issue is to assign time sequence in analysis, and the
effective stress from previously step will be used as reference to determine at each
location 1“ whether the current step is loading or unloading based on the integrated
stress value. Note the initial stress field is known as we require p,- = p0 = 0 so as to
produce a trivial stress field, and the stress solution can be computed step by step.
The next section will use numerical iteration to find the approximation of the solution

to equation (5.33).

108

5.3.2 Iterative Finite Difference Method

Since equation (5.33) is nonlinear, no analytical solution is available. We discretize
this equation and solve for an approximation to the exact transformation strains
by numerical iteration. Different iteration procedures can be constructed, and here
we will use the multi-variable Newton’s iteration. By meshing the field into discrete
points, the field variables are approximated by the field values at those points, and the
unknown variables are transformation strains at these points. The iteration procedure

is then used to solve a nonlinear algebra equation system.

Discretization

For the problem under consideration here, the simplest way of discretization is to place
equally-spaced points between 1‘,- and 7‘0. Let the number of total points (including
1* = 1",- and 1" = m) be N(N > 2), then the step size is h = (r0 —- ri)/(N — 1). The

coordinates of those N points are

rn=ri+(n—1)h, (n=1,2,...,N). (5.34)

For a field variable 1M1“), we use 11m to denote its value at m, i. e., w” = 1b(1'n)
(n = 1, 2, . . . , N) If the symbol of a field variable has subscripts, the symbol is put in
parentheses to avoid confusion, for example, eff at 1~ = Tn will be denoted as (53%;.
For the convenience of computer programming, we arrange the values of a field vari-
able at the discrete points in vector form. As a convention, bold symbols represent

vectors and matrices.

109

Radial transformation strains at those mesh points are

e’" = {(5% <e$T)2.....<e$’")N}T. (5.35)

and hoop transformation strains at those mesh points are

56 = {(53r)11(52r)2w ' '1(€t0r)N}T' (5'36)

Radial stresses at those mesh points are

a’" = {(a..)1.<a..>2,...,(0mm? (5.37)

and hoop stresses at those mesh points are

06 = {(099)1, (099)2A - - . . (GoalNlT- (5-38)

The integration H (r) and I (r) in equation (5.23) and (5.27) at rn can be approxi—

mated using the mid-point rule as

 

Hn = H(7‘n) = g: WHEAT) + (Effljl + Tj—lKEfrlj—l + (SETH—ll}
j=2 (5.39)
+r%(eg’">1 flag)...
and
In = I(1"n) = Z 5%] + Hg-l , (5.40)

where n = 1, 2, . . . ,N. Integral vectors H and I are then
H={H1,H2,...,HN}T (5.41)

and

I = {11,12,...,1N}T. (5.42)

With the above approximation of integration H (1") and I (1"), the constants Cg and C1

can be found as 2 2
_ W1 - (Po + 1Ner

 

 

 

2 2 ’
r — 1‘
N 1 2 2 (5.43)
CI ___ (P0 — 191+ INl'rer.
Ag, _ A;
At 1" = rn, the radial and hoop stresses are
(awn = Co + c, M, + In,
(5.44)
(099).. = co + (HA — cA>/r?. + I...
and the radial and hoop transformation strains are
2 _
(55% = (”0" (099)" aén.
2011,
(5.45)

 

2 n— rrn
(eight: (006)20n(0 ) 0:an

wheren=1,2,...,N.

111

If Pi and p0 are given and the stresses from the previous load step are given, the
discretized field variables are functions of 5" and 86 and the dependence on 1" is implied
through index n. We denote the functions for retrieving the above transformation

strains as

(n =1,2,...,N) (5.46)

Iteration Procedure

Suppose we assign an initial transformation strain field (ET)0 and (5:6)0 (the choosing
of initial value will be discussed in the sequel), an iteration scheme can be constructed

as

(1' =0,1,...) (5.47)

where superscripts denote the iteration number. Although this procedure is simple
and easy to implement, in practice the iteration procedure is prone to diverge except
for very good (15")0 and (€0)0. The Newton’s iteration with linear searches and
backtracking (Press at al., 2002) turns out to be very robust in solving this discretized
equation, and it usually produces converging results even for a initial transformation
strain field quite different from the converged solution.

The Newton’s iteration is widely used for root finding for nonlinear equations. It has

a quadratic convergence rate when the staring point for iteration is close to the real

112

solution. Suppose we want to find a root for a smooth nonlinear equation
F(:r) = 0. (5.48)

Let :1: = 1:0 be a first approximation of the real solution, the iteration procedure can

be constructed as
52"” = xi — F(:Ai)/F’(zi) (1': 1,2,. . .). (5.49)

Usually the iteration repeats until Irri‘l'l -- xil < ea; or IF (x)| < ef, where em and 8f
are two constants chosen depending on the given equation and the desired accuracy
goal. For a multi—dimensional situation where we want to find a root for the following
equation

F(a:) = 0, (5.50)

mi+1=mi— [1(a)] F(:ri) (1:1,2,...), (5.51)

0

where the initial point a: is given and the Jacobian J (m2) is

J(a:i) = g(mi). (5.52)

113

One commonly used criteria for stopping the iteration is that the difference between

two approximations is small enough
lmi‘l’l — xilg < ea; (5.53)

or the function is close to zero

|F(a:)|2 < ef, (5.54)

where l2-norm is used as the scalar measure of a vector.

However, the quadratic convergence of iteration (5.49) and (5.51) is a local property
of the given equation, and the success of these iterations is contingent to the selection
of initial value. Linear searches and backtracking (Press et al., 2002) can be use
to improve the convergence property. The basic idea is that each iteration should
“improve” the m, and the criteria chosen for determining a better a: is to decrease the

value of

f(a:) = F(a:) -F(a:)/2. (5.55)

Define the Newton step

. . -1 .
5.5% = — [Jag] Fat), (5.56)

then the modified Newton’s iteration will be

xi“ = xi + A613 (0 < A g 1). (5.57)

114

For each iteration, first let /\ = 1 and perform a normal Newton’s step. Then we

check if f(a:i+1) < f ((1:2), if this is true, then we proceed next iteration; otherwise,

the value of A will be reduced until f (2:241) < f ($2) is satisfied. Note we have

Vf-6a:=(F.J)o(—J—1-F)=—F-F<O,

which guarantees such a A can be found.
Let

a: = as")? (59107"

and

E(:.:) = {E{(a;), E§(a:), . . . ,E;V(m),E§(m),E§(a.-), . . . ,E?V(m)}T.

Equation (5.46) can be written as

F(a:) = a: - E(a:) = 0.

Starting from a initial (130, we need to find :1:* such that

F(a:*) = 0.

(5.58)

(5.59)

(5.60)

(5.61)

(5.62)

For each time step, we will use the resulting a: from the previous time step as the

initial value for iteration. For the first time step, we use 1:0

115

= O as the initial value

for iteration. Note for this specific form of F(:z:), F(a:z) = mi — E(:z:z) = sci — (CH-1,
therefore if we chose em = ef for equation (5.53) and (5.54), these two convergence

criteria are in fact identical.

Computation of J acobian

From equation (5.61), the Jacobian for this equation can be computed as

_ 8F _ 1 0E(a:)

J(a:)—5;— 0:1: '

 

(5.63)

The following gives the expressions for computing this J acobian. As a convention, we
will use two subscripts to represent the element of matrix, for example, Mi, j denotes
the 2th row and jth column element of matrix M.

Define H ’7' = BH/Ber, H ’0 = BH/Beg, both H 1" and H16 are N x N matrices.

From equation (5.39), the components of H17" and H19 can be computed as

10 fori=10rj>i,

Eh
—1‘1 forj=1and17£1,

9r _
HM — ( Eh (5.64)
forj=iand17é1,

 

L Ehrj for1<j<1',

and

116

0 fori=1orj>i,

Eh

-—2—1*1+E1‘% forj=1and17é1,
H130]. =

’ Eh

._ 2 .:. .
-—2—1"J ET] for} zandzaél,

(5.65)

 

Ehrj for 1 <j <1',
where 1',j =1,2...,N.
Define I ’7‘ = BI/aer, 1’9 = 81/659, both I 1” and I10 are N x N matrices. From

equation (5.40), the components of I ’7' and 1’6 are

 

 

' 1' 1'
h 1 Hi . Hf ._1
1T _ _ 1.7 1.7
1233' _ 2 2 7.3 "' 7.3 1 (5-66)
l=max(j,2) l l—1
and 0 6
i H’ . H’ .
16 _ h 1,] 1,3—1
1m — 5 Z ,.—3 + 7.3 ’ (5‘67)
l=maas(j,2) l 1‘1

where 1',j= 1,2,...,N.

117

Define 071” = 007/057”, 0"1 =0ar/0e, 06 000 /01-: 0'6 =0o'0/05, all of

which are N x N matrices. From equation (5.44), they are calculated by

 

 

 

 

 

r
T17: (TZ2— T2)TI2VT 117‘ + 11:7
2.1 (for ‘ If”? N’j
T16 __ (T22 — 7%)7‘12V 116 + 110
N 1 z
< 2 2 2 (5.68)
061T (Ti + r1)TNTI1T +T_H1:71+ I1j
2.] (va _ .37.? NJ ,
2 2 2
61.0 __ (7'2 + T1)TN I16 +_;H19j + 116
l 219 (1‘12V — 7%)112 N1—7+ 1'2.
where 1’, j = 1,2, . . . , N. Note in equation (5.68), 1*,- stands for the radial coordinate

for 1th discrete point rather than the inner radius.

Note H ,1, H10, I17” and I10 are all lower triangular matrices. More importantly,
H17, H10,I1", I10, 0’3", 0’39, 001’", 0016 are all the same for all iterations for a given
mesh, therefore they are required to compute only once if we fix the mesh (usually

the case).

118

Define 6"” = 0ET/0e’", 5’36 = BET/366, 59”” = 0E9/0e’", 50,9 = 0E9/059, all of

which are N x N matrices. Their components are calculated by

137' AU 7‘,r 300,1'
zyjzxiaz Ji+X ”2'0

A0130 30 0,0
52731.7 =X20’z +Xz 07:17
J (5.69)

60,1"

Eij 2X15 07‘,T+X1:C'0_6,T

i,j

 

0,0 B0130 000,0
K 67:17.: —Xi0 +Xi 21.7.

where

r ._ a .2 .
X54 _ $961+ [2(0rr72 ( 00);] (€/_ i”,

 

202' 0’2'
[7 (07‘7‘7i(000)i ‘401'2l I 51‘
< X2._ =2_:{— g,+ 2% (5,—07)}, (5.70)

 

 

0 ""0101“2 I i
Xc_ Tim” [2(00)z ( )l «VIE?»

20,-
and i,j = 1,2,...,N. Note here 5,; = g(a,,(gs),-) and (g = g,a(a,-,(gs),-). When
evaluating equation (5.69), the effective stress (am-)- and (099),- has been computed
by (5.44). Therefore the effective stress 02' can be inferred and the loading type can
be obtained by simplify comparing this 0,; to the o,- of previous load step. The switch
value (gs), can be computed based on either (2.9) or (2.11) depending on the loading
type and the g0 and so are martensitic volume fraction and effective stress at 1",; from

previous load step. The value of .5,- and {15 can then be easily obtained.

119

Finally, the Jacobian J can be written in the following block matrix form

J = 1 — . (5.71)

For no transformation strain situation, Q = E;- = 0 hence X 2’4 = X 2B = X 20 = 0, and

immediately J = l.

5.4 Numerical Results

A program implementing the iteration algorithm described in Section 5.3 is developed
using Mathematica, a scientific computing software (Wolfram, 2003). We denote this
program as iterating program. The responses of the plates under arbitrary traction
history can be solved by the iterating program. To verify the accuracy of this program,
the numerical solutions obtained by the iterating program, which will be referred to as
iterating solutions, are compared with the solutions obtained by the shooting method
described in section 5.2, which will be referred to as shooting solutions. First, we
show the the convergence of the iterating solution upon increasing the mesh points.
Then we simulate a plate subject to several different traction histories and present
the numerical results including the stress and martensitic fraction distributions, the
stress history and martensitic traction evolution at different locations, the movement

of phase interfaces, and so on.

120

 

5.4.1 Comparisons with the Shooting Solution

Here we use the iterating program to solve for the responses of the plate under the
same traction history and time step as those in Section 5.2. Although from the
shooting solutions the plate subject to the traction history will have global loading
responses followed by global unloading responses, such information is not required
for obtaining the iterating solution as the program will update the loading type au-
tomatically.

The shooting program is also implemented in Mathematica. Some implementation
details on the shooting program and the iterating program are provided here before
we compare the results from these two programs. The shooting program involves
numerical integration of a given ODE and root finding. Numerical integration is
given by NDSolve of Mathematica with accuracy options set to AccuracyGoal -—1 8
and PrecisionGoaI —> 8; and root finding is given by FindRoot of Mathematica with
accuracy options set to AccuracyGoal —> 8 and PrecisionGoal —+ 8. For options
AccuracyGoal —1 a and PrecisionGoal —> p, Mathematica attempts to maintain the
numerical error for a given quantity :1: to be less than 10‘“ + lelO‘p (for details,
please refer to Wolfram (2003)). For the iterating program, we use the convergence
criteria ex = ef = 1.0 x 10‘16, where ex and 6f are defined in equation (5.53) and
(5.54) respectively. This will require the solution found to satisfy |F(m*)|oo < 10—8,
which should be accurate enough for our purpose. Note for the iterating solution,
the field variables are given only at mesh points, while for the shooting solution,

NDSolve gives an interpolating function which can give all the values between r,- and

121

r0 numerically. The comparison will then be made only at the mesh points used in
the iterating program.

Figure 5.4 gives the stress distributions for the traction history (5.15) at t = 1 (Pi =
500 and p0 = 300). From t = 0 to t = 1, the plate is globally loaded to an ACM phase
distribution. The solid curves correspond to the stresses from the shooting solutions,
and the scatter points correspond to the iterating solutions. Note even with mesh
points of only 5, 9, and 17 points, the iterating solutions are already very close to the
reference solution obtained by the shooting method. With a mesh of 5 points, the
numerical results at those points (solid points in the figure) differ slightly with the
shooting solution. For more mesh points, the differences are hardly noticeable in the
plot.

Figure 5.5 gives the stress distribution for the traction history (5.15) at t = 1.5,
hence Pi = 2.500 and p0 = 1.500. From t = 1 to t = 1.5, the plate undergoes global
unloading. It keeps an ACM phase distribution although reverse transformation has
taken place in the plate. The solid curves correspond to the stresses from the shooting
solutions, and the scatter points correspond to the stresses from the iterating solution.
Similarly, the iterating solutions are very close to the reference solution obtained by
the shooting method and a larger number of mesh points produces more accurate
solutions.

To determine the mesh points needed for simulation, we need to estimate the errors

of the iterating solution. At a given time, let the field value from iterating program

122

 
 
   
 

 

 

' Spoints
C] 9points
4 ‘ Cl 17 points
2 _
10 12 1.4 1.6 18 20
e 1 l 1 J 1 I J l 1 I r
W
0'
_2 _ 00
-4 — Orr
-6 _

Figure 5.4: Comparison of stress distributions (normalized by 00) between the it-
erating solution and the shooting solution at t = 1. Here E = 50 GPa, a = 0.05,

of” = 00 = 100MPa, 0§T = 200, GET = 300, of?" = 400, and r0 = 2r,. The

applied tractions correspond to t = 1 are p,- = 500 and p0 = 300. The solid curves
correspond to the shooting solution, while the scatter points are results from the
iterating solution with mesh points of 5, 9, and 17 respectively.

be 052' and the result from shooting be ¢(r). The absolute error is defined as
errabs = max(|d>,- — ¢(ri)|), for i = 1, 2, . . . , N. (5.72)
The relative error is defined by

errrel = errabs/(z—fi, (5.73)

123

. I 5 points
‘ CI 9 points
El 17 points

 
   
 

 

 

 

Figure 5.5: Comparison of stress distributions (normalized by 00) between the iter-
atigg solution and the the shooting solution at t = 1.5. Here E = 50 GPa, a = 0.05,
of = 00 = 100MPa, 0.5271 = 200, JET = 300, of? = 400, and r0 = 2r,-. The
applied tractions correspond to t = 1.5 are pi = 2.500 and p0 = 1.500. The solid
curves correspond to the shooting solution, while the scatter points are results from
the iterating solution with mesh points of 5, 9, and 17 points respectively.

where 03 is a reference magnitude of function (0(r) defined by

 

fl," |¢3(0)Id0

7'0—7'7:

03 = (5.74)
Table 5.1 lists the absolute errors and relative errors for stresses at t = 1, and Table

5.2 lists the absolute errors and relative errors for stresses at t = 1.5. These two

tables demonstrate that the iterating solutions converge to the shooting solution as

124

 

the number of mesh points increases. The relative errors for 009 are larger than
those on». This may be caused by the following two reasons: (1) Evaluating 066
involves both Hn and In, but evaluating or,» only involves In. Since both Hn and
In are integrations approximated by the mid-point scheme, more truncation errors
can be introduced in evaluating 099. (2) Constants co and c1 in equation (5.43) are
solved by matching boundary conditions for or,- instead of 099, therefore the errors
for arr at r, and r0 are zero. This can be an advantage for arr in the context of
error estimation. On the other hand, all the errors are quite small and the difference
between the iterating solutions and the shooting solutions is hardly noticeable in
Figure 5.4 and 5.5. The errors decrease super-linearly upon increasing the number of
the mesh points, and they are negligible (less than 0.01%) when 129 mesh points are
used.

The responses of the plate are solved for quite a few other traction histories that
produce global loading and global unloading solution using the iterating program and
the shooting program. All those traction histories are proportional loading paths, and
have the following expression:

152-(t) = Hitoo and 150(t) = Hotao for 0 S t S 1,
(5.75)

13,-(t) = 112-(2 — t)00 and 150(t) = H0(2 — t)00 for 1 < t S 2,

where H,- and H0 are two constants for a specific traction history and again 00 = 0;”.

These paths produce global loading solutions for 0 S t S 1, and global unloading

solutions for 1 g t S 2. Traction histories with {11,3110} as {0,2}, {2,4}, {3,4},

125

 

Mesh Point Number 5 9 17 33 129

 

errabs for arr 1.5E — 2 3.6E — 3 9.5E — 4 2.4E — 4 1.5E — 5
errabs for 066 9.6E — 2 2.6E — 2 6.1E — 3 1.3E — 3 8.8E — 5
errabs for 0 4.2E — 2 1.1E — 2 2.7E —- 3 5.8E — 4 4.0E — 5
errrel for arr 0.391% 0.095% 0.025% 0.006% 0.0004%
errrel for 099 9.579% 2.577% 0.608% 0.130% 0.0088%
errrel for a 1.240% 0.338% 0.080% 0.017% 0.0012%

 

Table 5.1: Absolute errors and relative errors at t = 1 for the iterating solutions
compared to the shooting solutions.

 

 

Mesh Point Number 5 9 17 33 129

errabs for 07-7- 9.5E — 3 2.4E - 3 6.5E — 4 1.6E — 4 1.0E — 5
errabs for 099 4.7E — 2 1.3E — 2 3.0E — 3 6.1E — 4 4.2E — 5
errabs for a 1.9E — 2 5.0E — 3 1.2E — 3 2.4E — 4 1.7E — 5
errrel for arr 0.495% 0.127% 0.034% 0.008% 0.0005%
errrel for 099 9.474% 2.537% 0.594% 0.123% 0.0085%
errrel for a 1.076% 0.293% 0.069% 0.014% 0.0010%

 

Table 5.2: Absolute errors and relative errors at t = 1.5 for the iterating solutions
compared to the shooting solutions.

{6,3}, {8, 2}, {8,5}, {10,0}, et al., are solved. The comparison between the results
from these two programs again shows that the iterating solutions are accurate and
efficient. Since the results are similar, they are omitted here for conciseness.

Table 5.3 shows the estimated errors for traction history (5.15) with mesh points
of 101 (the step size in r direction is h = O.01(r0 — r,)) for all time steps. The
maximum relative error for all these time steps is around 0.0146%, which is very
small. Usually a number of 101 mesh points is sufficient for obtaining a quite accurate
stress distribution, however, to obtain the movement of phase interface may require
more refined mesh, which will be discussed later.

The number of time steps used is also very important in finding both the shooting
solutions and the iterating solutions. For the global loading solution and global

unloading solution, however, it does not affect the results at a given time as long as

126

the time at which the loading type changes is included. For instance, in this example
the plate changes from global loading to global unloading at t = 1, therefore t = 1
must be included in time steps. For generic traction histories, the selection of time
step can be of great impact on the solutions. A practical approach is to use crude time
steps first, and refine the time step until the solutions at the desired time converge.
Although an automatic scheme to adjust time step is possible, we will choose the

time step manually in this study.

5.4.2 Numerical Results for Plates Subject to Arbitrary Trac-
tion History

This section will show the iterating solutions for a few non-proportional traction
histories where the shooting program is not applicable. In what follows, we will only
study an annular plate with r7; = 2r0, and the material parameters are E = 50 GPa,
a = 0.05, 0;” = 00 = 100MPa, 05F = 200, 0571 = 300, of? = 400, and r0 = 2r,.
Simulations with different material constants and plate geometries can be performed

easily but they do not reveal much qualitative difference.

127

Example 1

This example shows the responses of the plate under the load history

7
13,-(t) = 2t00 and 150(t) = 2t00 for 0 S t S 1,

131'“) = 200 and 15005) = (1 +000 for 1< t S 2,
( (5.76)

13,-(t) = too and 150(t) = (5 — t)00 for 2 < t S 3,

 

0i“) = (6 - t)00 and [30(t) = 200 for 3 < t S 4.

k

This traction history will be referred to as path I in the sequel, and the traction
history is plotted in Figure 5.6 together with part of the loading structure map. In
(19,3190) plane, the traction follows path 0 —+ a -—+ b —+ c —+ a. Note at t = 1 and t = 4,
we have p,- = p0 = 200. Thin curves in the figure are part of the loading structure
map, and for path 0 —> a, the plate remains .A phase and an equi-biaxial stress field;
for what follows, we can not determine the phase type only by the structure maps
since reverse transformation may be involved.

Figure 5.7 is the plot for effective stresses at t = 1 and t = 4, and Figure 5.8 is the
plot for radial stresses and hoop stresses at t = 1 and t = 4. As mentioned above, at
t = 1, a = 200, which corresponds to the dashed line in the figure. The scatter points
and the solid line are effective stress at t = 4. The solid dot points are computed with
a total of 40 time steps (At = 0.1) and 33 mesh points in r direction; the hollow dot
points are computed with a total of 80 time steps (At = 0.05) and 33 mesh points in r
direction; the solid curves are computed with a total of 80 time steps (At = 0.05) and

101 mesh points in r direction, and the 101 discrete points are connected together.

128

 

)po/O'jfl
a— (2, 2)
5 * b— (2. 3) *
c— (3, 2)
4 _ _

 

. i a . .pi/Gfr
1 2 3 4 5 6

 

Figure 5.6: Traction history of of Example 1. The load follows a path 0 —+ a —+ b ——>
c —> a in (11,3190) plane. Other thin curves are part of the loading structure map. This
load history will be referred to as path I.

From these two figures, we find that doubling the number of time steps and increasing
the mesh points has relative little effect on the stresses at t = 4, therefore, we will use
a step number of 80 and 101 mesh points to compute the solutions. In what follows,
we assume that the number of time steps and the number of mesh points chosen are
sufficient to satisfy our accuracy requirements, therefore no further comparisons on
solutions corresponding to different numbers of time steps and mesh points will be
given.

At t = 1, the plate will have p,- = p0 = 200 and a uniform stress field as arr(r) =
099(r) = —200 and 0(r) = 200. Although at t = 4 we also have p; = p0 = 200 but
the stress field is not uniform, as shown in Figure 5.7 and 5.8. Starting at t = 1,
after a load cycle following a -) b —-> c —> a on (pi, p0) plane (Figure 5.6), the stresses

do not return to the uniform stresses at t = 1. Examining the martensitic fraction

129

I t=4 (40 steps, 33 mesh points)
D t=4 (80 steps, 33 mesh points)
11:4 (80 steps, 101 mesh points)

    

 

 

1.8 T I I I l I T I I r 7‘

Figure 5.7: Comparison of effective stress distributions (normalized by 00) between
t = 1 and t = 4 for to traction path I. The applied tractions at t = 1 and t = 4 are
both p,- = p0 = 200. The scatter points and solid curve are effective stresses at t = 4.
The dashed-line is the effective stress at t = 1 and 0(r) = 200.

reveals that the whole plate has been transformed from originally A phase to (‘3 phase.
This is an example showing the history dependency in plate responses.

Figure 5.9 shows the effective stress history for r = 1, r = 1.5 and r = 2. For
0 S t g 1, the plate has a uniform effective stress field and the effective stress keeps
increasing. For 1 < t g 2, forward transformation takes place, and everywhere in the
plate is loading. For 2 < t S 3, initially all these three locations are unloading, then
the outer part of the plate keeps unloading, while the inner part of the plate starts
loading. No forward transformation takes place at r = 1 since the effective stress is
below JET. At r = 1 and r = 1.5, the reverse transformation is activated when t
approaches t = 3. For 3 < t S 4, the inner region of the plate switches from loading

to unloading, and the outer region of the plate switches from unloading to loading.

Locations in the middle of the plate at first continues to unload, but then it switches

130

 

I 11:4 (40 steps, 33 mesh points)
“ El t=4 (80 steps, 33 mesh points)
t=4 (80 steps, 101 mesh points)

    

 

 

Figure 5.8: Comparison of radial stress and hoop stress distributions (normalized by
00) between t = 1 and t = 4 for traction path I. The applied tractions at t = 1 and
t = 4 are the same, with p,- = p0 = 200. The scatter points and solid curves are
stresses at t = 4. The dashed—line is the stress at t = 1 and an (r) = 099 (r) = —200.

to loading. During 2 S t S 4, we observe that the effective stress does not have the
property of monotonously decreasing along r direction and loading and unloading can
occur at the same time. Figure 5.10 shows the £ — a plot for that r = 1, r = 1.5 and
r = 2 for this traction history. Note the evolution of martensitic fractions for these
three locations follow the model described in Chapter 2.

Specifically for 2 g t S 3, from Figure 5.9 we find that there is loading and unloading
coexisting on the plate when t approaches t = 3. Figure 5.11 shows the effective
stress distributions for the time steps with 2 g t S 3. In the figure, the tOp most
curve corresponds to the effective stress at t = 2, and the curves on the lower part
are those with t close to 3. The difference in effective stresses between two consec-
utive steps (Ad[t] = o[t] — o[t — At]) are plotted in Figure 5.12. It is shown in the

figure that everywhere in the plate is unloading for 2 S t < 2.7, but there is no

131

 

 

 

 

4 — o

- ---r=L5

T —— 1':
3 2% ..................................................
\

‘ 0Jar
2 -‘ ----- '5' """""""""""""""""" J """""""""""" ‘V ‘c‘;; """"
1 _
O 1 j I i I T I 1 t

0 1 2 3 4

Figure 5.9: Time history plot of effective stresses (normalized by 00) at r = 1, r = 1.5
and r = 2. Applied traction is path I, and a time step At = 0.05 and 101 mesh points
are used in the computation.

reverse transformation activated as o > 0?? Similarly, at t = 2.7, t = 2.75, and
t = 2.8, although in part of the plate the effective stress continues to decrease, no
reverse transformation takes place as the effective stress is everywhere above 0?.
Figure 5.13 shows the difference in martensitic fraction between two consecutive steps
(A45 [t] = 5 [t] — E [t — At]) for 2 < t < 3. There is no change in martensitic fraction
until t = 2.85.

Similarly, for 3 S t _<_ 4, from Figure 5.14 we again find that there will be loading and
unloading coexisting on the plate. Figure 5.14 shows the effective stress distribution
for 3 S t S 4. Throughout this time period, the inner part of the the plate undergoes
unloading, while the outer part of the plate undergoes loading. Figure 5.15 shows the
difference in effective stresses between two consecutive steps. Recall that to initiate

forward transformation requires a 2 GET and to initiate reverse transformation

requires a g 0,571. At t = 3.05, there will be reverse transformation taking place

132

   
 

 

 

 

1 o . _ x _,.=, . ----------------------------------
I! ”
—o—1‘=l.5 '1' '1’
__ __ ... ' 1
D 7' _2 1' WX-x-X—x—
1"
0 5 ‘1 moooooo-oooooo
—( .
. I :OODDOOOOUCODOCCOD
I"
0
1
I
a
I
I:
I, 0/60
0.0 A22:=:::::::q:.‘:.:::::::::::::T I I
0 l 2 3 4

Figure 5.10: Evolution of martensitic fractions at r = 1, r = 1.5 and r = 2. Effective
stress a is normalized by 00. Applied traction is path I, and a time step At = 0.05

and 101 mesh points are used in computation.

for ra < r < rb (please refer to Figure 5.14 and 5.15); when r < ra the effective

stress is decreasing but a > 051‘; when r > rb the effective stress is increasing but

a < GET. As a result, at t = 3.05 reverse transformation takes place for ra < r <b,
while for other locations, neither forward transformation nor reverse transformation
takes place. Note for locations with r < ra and r > rb, although there is no change in
martensitic volume fraction for this time step, the transformation strain can change
due to the change in stresses. Figure 5.16 shows the difference in martensitic fraction

between consecutive steps (A{ [t] = g [t] — 5 [t — At]). It is shown in the figure that

the section of the plate undergoing reverse transformation moves inwards as time

increases.

Example 2

This example uses a traction history that linearly changes between two points on
{pl-,po} plane. During the first second, the inner and outer tractions {19,3190} increase

linearly from {0,0} to {400, 200}; after that, the tractions oscillate linearly between

133

 

3.5 —

.---—------------------------------------—------.

 

---
'up
”-

3.0 """ (I: _ 0’?

 

--—
-’-‘--
-
‘-—-
‘
” ----
"‘ ------------—--
-
" --------
av-
‘-
—-
""
’
‘ ---
“‘ ----------‘---
--------
----
---
---
“
“'-
------‘----‘----‘
-------‘-----
----------

2.5 "‘: V‘—
‘&

 

 

 

 

 

 

 

‘ W“ ------------------------------------------
----::::: -------------- - -----r T
2.0 \W .,-- ---;og
- i=3 t=2.95 """""
t=2 90 t=2.85 t=2. 75
15 t—2.80 ’r
o I I I l l I 1 l 1 I
1.0 1.2 1.4 1.6 1.8 2.0

Figure 5.11: Effective stresses (normalized by 00) for 2 S t S 3. Applied traction
is path I, and a time step At = 0.05 and 101 mesh points are used in computation.
During this period, no forward transformation takes place. Part of the plate undergoes
reverse transformation for 2.85 S t S 3.

{400, 200} and {200, 400} for a period of two seconds, as shown in Figure 5.17. This
traction history will be referred to as path 11 in the sequel. In the computation, we
calculate the time period 0 S t S 7. The time step is set to At = 0.01, and 101 mesh
points are used in r direction.

Figure 5.18 shows the changes of effective stresses about time for three locations
on the plate, r = 1, r = 1.5, and r = 2. It is found that when t Z 2 we have
0(t + 2k) = 0(t). This is due to the fact that at time t = 2k (k = 1,2,. . . ), the whole
plate is transformed into pure martensite. As mentioned in the discussions on 9,; in
Section 5.1, previous history does not affect the future responses since gs(r) E 1.
Then such a pattern repeats and the ensuing responses are periodic functions of time

with a period of two seconds. Note when t approaches t = 3, the outer part of the

134

 

 

 

 

 

 

1.0 1.2 1.4 1.6 1.8 2.0

Figure 5.12: Difference in effective stresses (normalized by 00) between two consec-
utive steps (A0[t] = a[t] — 0[t — At]) for the period 2 < t S 3. Applied traction is
path I, and a time step At = 0.05 and 101 mesh points are used in computation.

plate has effective stress below 0,511 and is unloading, therefore reverse transformation

takes place.

Example 3

This example again uses another traction history that linearly changes between two
points on {p33 p0} plane. During the first second, the inner and outer tractions {p1, 190}
increase linearly from {0,0} to {400, 200}; after that, the tractions oscillate linearly
between {400, 200} and {200, 300} for a period of two seconds, as shown in Figure
5.19. This traction history will be referred to as path III in the sequel. In the
computation, we calculate the time period from t = 0 to t = 20. The time step is set

to At = 0.01, and 201 mesh points are used in r direction.

135

 

 

 

 

0.01 —
2 < t <2.85
0.00
1
-0.01 4
-0.02 —
-0.03 —
'0.04 r I I I I I r I 1 IT
1.0 1.2 1.4 1.6 1.8 2.0

Figure 5.13: Difference in martensitic fractions between consecutive steps (A5 [t] =
E [t] — €[t — At]) for the period 2 < t S 3. Applied traction is path I, and a time step
At = 0.05 and 101 mesh points are used in computation.

Figure 5.20 is the effective stress history plot for r = 1.0, 1.2, . . . ,2.0. Initially effective
stresses increase for 0 S t S 1, then the effective stresses oscillate as time proceeds.
The effective stress responses of the plate in a period of 2 seconds gradually stabilize
as the load cycle continues. Figure 5.21 shows the evolution of martensitic fractions
for r = 1.0,1.2,.. . ,2.0. These stress history plots demonstrate the stabilization of
the stress responses on the plate. The plots are quite similar to Figure 2.6, in which
the effective stress changes between two fixed values. However, here the effective
stresses do not change between two fixed values, instead, the effective stresses turn
to change between two converging values as load cycle continues.

Figure 5.22(a) singles out effective stresses at t = 1,3,...,19. Effective stress dis-
tributions at t = 1 and t = 3 are quite different, but as time proceeds, the effective
stress at t = 2k — 1 (k = 3,4, . . . ,9) converges rapidly. Figure 5.22(b) singles out
effective stresses at t = 2, 4,... ,20. Similarly, effective stress distributions at t = 2

and t = 4 are quite different, but as time proceeds, the effective stress at t = 2k

136

 

 

 

3.0 a?
k
‘, t=3
1 \\
\ \\
\
2.5 —( \\ \\
t \
\ \ \
\ \ \
\ s \
~ \\ \\ \\
q \“ \\ \‘ \\
\ \ \ \
p \ \ \
\ ‘\ s \ \ t_4
\‘ ‘ \\ \\ \\ _
~~ \‘ \‘\ \\ \\ \\
~ ‘s‘ \\“‘\“\\:\‘\\‘ 7.0—142 .......... -
2.0 ... ‘ \‘ ‘ ~ ‘ ‘s ‘\ """""""" 0
‘~\ $3:- W----_—: s
..-_ _-_ - —- ~ - ,r.~? :‘3‘~. ---_"_’_‘ ..................
.1 -
7.
1'5 I I I I l I I I I I
1.0 1.2 1.4 1.6 1.8 2.0

Figure 5.14: Effective stresses (normalized by 00) for the period 3 S t S 4. Applied
traction is path I, and a time step At = 0.05 and 101 mesh points are used in
computation.

(k = 3, 4, . . . , 10) converges rapidly. Although this figure only gives the plot for ef-
fective stress, other stress components and the martensitic fraction exhibit similar

converging trends.

Example 4

This example shows the movement of phase interfaces. The load history is expressed

as

132-(t) = 4t00 and 150(t) = 0 for 0 S t S 1,
0i“) = 4(2 “ ””0 and 750(t) = 4(t — 1)0’0 for 1 < t S 2, (5-77)

Fifi) = 0 and 730“) = 4(3 - t)00 for 2 < t S 3.

137

 

- 0.05 —
' XII/1’ Unloading

 

 

‘0. 10 I I I I I I I I I Ir
1.0 1. 2 1.4 1.6 l. 8 2.0

Figure 5.15: Difference in effective stresses (normalized by 00) between consecutive
steps (Aa[t] = a[t] — a[t — At]) for the period 3 < t S 4. Applied traction is path I,
and a time step At = 0.05 and 101 mesh points are used in computation.

This traction history will be referred to as path I V in the sequel. In the computation,
the time step is set to At = 0.01, and 201 mesh points are used in r direction.

In Chapter 3, we define the interface r A as the boundary between an inner (‘3 ring
and an outer A ring, and the interface TM as the boundary between an inner M ring
and an outer 63 ring. For the arbitrary traction histories, the non-increasing of the
effective stress along radial direction does not hold any more. The martensitic fraction
is not necessarily non-increasing along radial direction. However, in the numerical
simulations performed so far, the martensitic fraction is always non-increasing along
radial direction. Hence we again use r A to denote (“I/A interface and TM to denote
M/C‘f. Also in all our simulations, at most only one C/A interface and one M/C

interface are found.

138

 

ra=1.42 rb=1.75

 

 

 

 

'0.m6 I I I I F I I I I Ir
1.0 1.2 1.4 1.6 1.8 2.0

Figure 5.16: Difference in martensitic fractions between consecutive steps (Af [t] =
g [t] — 5 [t — At]) for the period 3 < t S 4. Applied traction is path I, and a time step
At = 0.05 and 101 mesh points are used in computation.

For the discretized mesh here, the (‘3 /.A interface is retrieved by

7'A = (rn + rn+1)/2, (5.78)

where 5n > 0 and {n+1 = 0. When fin = 0 for all the mesh points, there is no 6/04
interface in the plate; for convenience of plotting, we will set r A = Ti- On the other
hand, when {n > 0 for all the mesh points, there is no (‘3 /./1 interface in the plate and
we will set r A = r0.

The M/C interface is retrieved by

TM = (Tn + 7"n+1l/21 (5.79)

139

 

 

 

5 — pi/Go
1 ----Po/oo
4""'" ""',.I""‘ ""}S"'" ""75"".
J l ’l\ I ’l\ I ’l\ I
' 1 ' 1 : , : 1 : , : \ :
I I \ \ \
3 7 i I i , i i 1 i i
I I, I \ I I I \ I I I \ I
1 :1 : \:z : \L’ : 1):
2a _____ z .......... t .............. ----
[I I I I I I I
_ ’l I I I I I l
1’ 1 l I I i i I
1- , : : : : : : :
_ I I I I I I I I
z i I i I l i i
O . 1 . i . i . i . , . i . , t
0 l 2 3 4 5 6 7

Figure 5.17: Traction history for Example 2. This traction history is denoted as path
II.

12—o

 

.1

10—

.5 """""""""""" A 'A

6—

 

 

1 I I

-l : :

4 o : 1 n :
-—I ..................................... -4. ..... a- ..... db .......

1 I 1

. 1 l i

d A} '

\ I 1

2 — O... . A, ...: ......... . .......... .‘.- .q.
' I" “ e ‘ I". ‘3

1 I 1

-l a'

O T r r I: ' i I f j I f t

 

 

 

 

 

 

 

 

 

 

Figure 5.18: Time history plot of effective stresses (normalized by 00) at r = 1,
r = 1.5, and r = 2. Applied traction is path II, and a time step At = 0.01 and 101
mesh points are used in computation.

where En = 1 and €n+1 < 1. When {n < 1 for all the mesh points, there is no M/C
interface and we will set TIM = Ti- On the other hand, when 5n = 1 for all the mesh
points, there is no M/(i‘ interface and we will set TM = r0.

Figure 5.23 shows the movement of r A (C/A interface) and TM (MI/(‘3 interface) for
this traction history. Initially the plate is pure austenite, hence r A = TM = ri. At
t = 0.33 we have the appearance of C/A interface where the plate initiates forward

transformation at the inner boundary. This C/A interface keeps moving outwards

140

 

 

 

 

4 ... .............................................. ‘ --.
-I I I I I I I I I I
I I I I I I I I I I
3-«~-I--.,-+1.-4-1-I-1-I-r-I-1c4-r In. I— -'r-.
I ’1 \ I I I I I '1‘ I ’I‘ I ’I I I‘ I ’1‘ I ’5 I I
- I I I I I I I
it” i “I ,’ : X I I’ : ‘\ I’l’ : \ \l’l’ : ‘\‘l’[’ i ‘I ,’ i ‘\\:’I' : ‘\\ :1” i “i ,’ i
2 -- -r--. -‘I’-- --’t--.--‘I--.--I-- -+- .--“I’-- .°-I-- "fun-If".
"I I I I I I I I I I I I I I I I I I I I
1 '1 I I I I I I I I I I I I I I I I I I I
’1 I I I I I I I I I I I I I I I I I I I
1_ , I I I I I I I l I I I I I I I I I I I I
I I I I I I I I I I
,' . I . I . I . I I I . I I I I I I I . I
I I I I I I I I I I I I I I I I I I I I
I I I I I I I I I
O , I ; 4 j I 41' I ; I I + j I j I I I I I it
0 2 4 6 8 10 12 14 16 18 20

Figure 5.19: Traction history for Example 3. This traction history is denoted as path
III.

as a result of global loading, and the interface reaches r A = r0 at t = 0.62. As the
load keeps increasing, M/G interface emerges from inner boundary at t = 0.67. This
interface keeps moving outwards and the M ring expands.

At t = 1, the applied tractions switch their direction, and the movement of M/G
interface stops due to unloading. However, the reverse transformation is not activated
since the stress is above of”. Therefore, the 7171/63 interface stays at TM = 1.245 until
t = 1.35, when the reverse transformation is activated. The M/C‘? interface is then
moves back toward Ti; at t = 1.44 the interface TM reaches r,- and the plate is in
E state. The M/C‘Z starts again at t = 1.68, and the plate is then loaded to full
martensite at t = 1.86. During this procedure, C/A interface never appears inside
the plate.

From 2 S t S 3, the plate is globally unloaded from martensite, and the M/C
interface reappears at t = 2.57 and moves toward 7'2" It disappears at t = 2.78. The
(“T/A interface reappears at t = 2.86 and moves toward 7'2“ It disappears at t = 2.91

whereupon the whole plate is converted back to A phase.

141

Example 5

This example demonstrates the capability of the program in handling complex trac-
tion histories. Here we chose a traction history with {19,3190} expressed as the following
functions of time:

132'“) = tsin(2t)00, (0 _<. t S 27f) (5.80)

130(t) = tcos(2t)00/2.

Figure 5.24 shows the traction history on {19,3120} plane together with part of the
loading structure maps. This traction history will be referred to as path V in the
sequel. A time step At = 0.01 and 201 mesh points are used in computation.

Figure 5.25 shows the effective stress history for r = 101.2, .. . ,2.0. The effec-
tive stresses oscillate as the tractions change. Figure 5.26 shows the evolution of
martensitic fraction for r = 1.0, 1.2, . . . , 2.0, where partial reverse transformation and
forward transformation take place during the time period.

Figure 5.27 shows the movement of 1' .A (63/11 interface) and TM (M/G interface) for
this traction loading history. Initially, no transformation starts as the tractions are
inside the FAL on the loading structure map, hence TA = TM = 7'2" We have the
appearance of C/A interface at t = 1.18 (p,- = 0.83100 and p0 = —0.41900). This
(3 ring expands to TA = 1.0075, and quickly it is transformed back to austenite at
t = 1.76 (p,- = —0.65000 and p0 = —0.81800). The 6/11 interface appears again at

t = 2.1 (p7; = —1.83000 and p0 = —0.51500), and it is pushed outward. At t = 2.39

142

 

(p,- = —2.38500 and p0 = 0.08100), the G/A interface reaches To and no A phase
exists on the plate. No further C/Jl interface appears in this simulation.

The forward transformation continues and at t = 2.46 (192- = —2.40700 and p0 =
0.25400) interface M/(‘B appears, and this interface propagates outwards until t = 2.64
(p,- = —2.22600 and p0 = 0.71000). The interface keeps at this location until t = 3.04
(p,- = —O.61300 and p0 = 1.48900), after which the interface moves inwards and
disappears at t = 3.21 (pi = 0.43800 ande = 1.59000). The M/C interface reappears
at t = 3.72 (p,- = 3.40600 and p0 = 0.74800), and keeps propagating outwards until
t = 4.28 (p,- = 3.25700 and p0 = —1.38800), after which the interface stays until
t = 5.38 (p,- = —5.23100 and p0 = —0.62800). The M/C interface keeps propagating
outwards until t = 6.16 (pi = -1.50200 and p0 = 2.98700), and for 6.16 < t S 6.28

the interface remains at TM = 1.7075.

5.5 Discussions

The iterating program is shown to be efficient in solving the responses of the plate
under various traction histories. For the shooting program, if the loading type is
misidentified at one location due to numerical error during the procedure of inte-
grating the ODE system, the following integrating steps can be incorrect because
the numerical error can be accumulated as integration proceeds. However, for the
iterating program, the contribution of a load type misidentification at a location is

diminished as a large number of other points have correct loading types. For example,

143

from equation (5.43) we can see that the effect of a load type misidentification due
to numerical error on the value co and c1 will be marginal.

The computation of phase interface by equation (5.78) and (5.79) will limit the accu-
racy of interface evaluation to an error within h/ 2 (h is the step size in 1" direction).
To improve the accuracy of interface, one way is to increase the mesh point, which will
increase the simulation time significantly; another way is to interpolate the effective
stress and martensitic fraction between mesh points, and then find the locations that
are (“f/A interface and/ or G/A interface.

The iteration procedure can be easily extended to other material models. The deduc-
tion of stress solution in section 5.3.1 does not involve any actual constitutive model.
The only requirement is the decomposition of strain into elastic part and inelastic
part, as shown in (5.6), which is a very general assumption. The actual material
model is then used to compute Jacobian for numerical iteration. Therefore, it is pos-
sible to implement a more complicated material model with appropriate modifications

on Jacobian.

144

 

Orr

errabs
006

0'

0M

errrel
000

0'

 

t= 0.05
t= 0.10
t= 0.15
t= 0.20
t: 0.25
t= 0.30
t= 0.35
t= 0.40
t= 0.45
t= 0.50
t: 0.55
t= 0.60
t: 0.65
t: 0.70
t= 0.75
It: 0.80
t= 0.85
t= 0.90
t= 0.95
t: 1.00
t= 1.05
t= 1.10
t: 1.15
t: 1.20
t= 1.25
t= 1.30
t= 1.35
t: 1.40
t= 1.45
t: 1.50
t= 1.55
t= 1.60

= 1.65

= 1.70
t= 1.75
t= 1.80
t= 1.85
t= 1.90
t = 1.95
t= 2.00

 

15E-8
23E—8
28E—8
3AE-8
50E-8
53E—8
59E—8
66E—8
llE-7
12E-7
13E—7
40E—5
57E—5
63E—5
47E—5
29E—5
18E—5
31E—5
27E—5
25E—5
24E-5
23E—5
22E—5
2JE—5
2JE—5
2JE-5
20E-5
20E—5
18E—5
IVE—5
15E—5
89E—6
49E—6
56E—6
15E-5
99E—6
28E—8
23E—8
15E—8
000E0

30E—8
43E—8
52E—8
67E-8
51E-8
42E—8
49E—8
57E—8
97E-8
11E-7
12E—7
6AE—5
83E—5
92E—5
11E—4
10E—4
89E—5
10E—4
12E—4
14E—4
14E—4
13E—4
12E—4
12E—4
11E-4
10E—4
95E—5
87E—5
79E—5
69E—5
58E—5
38E—5
35E—5
29E—5
32E—5
20E—5
52E-8
43E—8
30E—8
000E0

18E-8
3JE—8
4JE—8
5JE—8
GJE—8
58E-8
66E—8
74E—8
92E—8
10E—7
10E—7
74E—5
98E—5
10E—4
7AE—5
4JE—5
10E—5
BJE-5
44E—5
65E-5
6JE—5
58E—5
54E—5
5lE—5
4JE—5
44E—5
4JE—5
3JE-5
33E—5
28E—5
20E—5
52E-6
47E—6
80E-6
24E—5
18E—5
4JE—8
3JE—8
18E—8
000E0

 

0.00007)
0.00007)
0.00007)
0.00007)
0.00007)
0.00007)
0.00007)
0.00007)
0.00007)
0.00007)
0.000070
0.001870
0.002470
0.00247)
0.001770
0.001070
0.000670
0.000970
0.000770
0.00077)
0.00077)
0.00077)
0.000770
0.00077)
0.00077)
0.000870

. 0.00087)

0.00097)
0.00097)
0.00097)
0.00097)
0.00067)
0.000470
0.000570
0.001670
0.001470
0.000070
0.000070
0.000070
0.000070

0.00007)
0.00007)
0.00007)
0.000070
0.000070
0.000070
0.000070
0.000070
0.000070
0.000070
0.000070
0.010570
0.012870
0.013170
0.015170
0.01287)
0.010470
0.01137)
0.012270
0.014470
0.014570
0.01457)
0.01457)
0.01467)
0.01457)
0.014670
0.01467)
0.01457)
0.01447)
0.013870
0.012970
0.009670
0.009970
0.00967)
0.013070
0.00967)
0.000070
0.000070
0.000070
0.00007)

0.000070
0.00007)
0.000070
0.00007)
0.00007)
0.00007)
0.00007)
0.00007)
0.000070
0.00007)
0.00007)
0.00377)
0.00457)
0.00447)
0.00297)
0.00157)
0.00037)
0.00107)
0.00147)
0.00197)
0.00197)
0.00197)
0.00197)
0.00197)
0.00197)
0.00187)
0.00187)
0.00187)
0.00187)
0.00167)
0.00137)
0.000470
0.00047)
0.00087)
0.00287)
0.002770
0.00007)
0.00007)
0.00007)
0.00007)

 

 

Table 5.3: Absolute errors and relative errors for the iterating solutions compared
to the shooting solutions. Time step used is At = 0.05 and totally 101 mesh points
(step size h = (r0 — Ti) / 100) are used for simulation.

145

 

 

 

 

 

(a) r =1.0 (b) r=1.2

 

 

 

 

 

 

 

 

 

-o
4- .........................
t ' t
II I O I I I I I I I I I1
16 20 0 4 8 12 16 20

(d)7"=1.6

o -o
4- ------------------------- 41h ------------------------
2—------ 2—------

‘ t ‘ t
O ' lfilj l ' l 'l O ' l ' l ' l ' I ' I
0 4 8 12 16 20 O 4 8 12 16 20

(e) 7'=1.8 (t) r=2.0

Figure 5.20: Time history plot of effective stresses (normalized by 00) at r =
1.0,1.2,.. . ,2.0. Applied traction is path III, and a time step At = 0.01 and 201
mesh points are used in computation.

146

 

 

 

 

 

 

(b) 'r =1.2

(a) r=1.0

 

 

 

 

 

(d) T=1.6

(c) r=1.4

 

 

 

 

 

 

(D r=2.0

(e) T’=1.8

. , 2.0. Effective stress

a is normalized by 00. Applied traction is path 111, and a time step At = 0.01 and

Figure 5.21: Evolution of martensitic fractions at r = 1.0, 1.2, . .
201 mesh points are used in computation.

147

 

I I I I I I I I I I I I I I I I I I I IT
1.0 1.2 1.4 1.6 1.8 2.0 1.0 1.2 1.4 1.6 1.8 2.0
(a) (b)

Figure 5.22: Convergence of effective stresses (normalized by 00). Applied traction
is path III, and a time step At = 0.01 and 201 mesh points are used in computation.
(a) Effective stresses at t = 1, 3, . . . , 19. (b) Effective stresses at t = 2, 4, . . . ,20.

3.0

 

1

 

 

 

 

 

 

 

 

 

 

 

 

O T F r I ' I ' I ' 7.
1.0 1.2 1.4 1.6 1.8 2.0

 

Figure 5.23: Movement of TA (C/A interface) and TM (CM/C interface). Applied
traction is path V, and a time step At = 0.01 and 201 mesh points are used in
computation.

148

 

 

 

 

 

I——-————— ——--—‘—

_5 l 1 1 L

-6 —4 -2

 

0
PI/oo

Figure 5.24: Traction history for Example 4. This traction history is denoted as path
V.

149

 

 

 

 

 

 

 

 

 

 

 

 

 

-0- 0'
24 a ..................... - 161‘
16 - 12 — -------------------- -
""""""""""" 8 74-------------- --- --
81W - -- -- - 4; _________ -- -
O I I I I I I t 0 - I I I j I I t
0 2 4 6 0 2 4 6
(a) r=10 (b) r=12
)6 ' G
6 — --------------------- - _ ______________________
I 4
4 _ _______________________ ..
I _ ___________
2 I -------------------- 2 -
t t
0 ' I I I ' I 0 ' I I ' I
0 2 4 6 0 2 4 6
(c) r=14 (d) r=16
- 0' - 0'
4 ----------------------- 4 _ ........................
2 I. ---------- —--- - --- 2 _. ........... ---- --..
" t ‘ t
O ' I ' I ' I 0 I I r I
0 2 4 6 0 2 4 6
(e) T=1.8 (t) r=2.o

Figure 5.25: Time history plot of effective stresses (normalized by 00) at r =
10,12, . . . ,2.0. Applied traction is path V, and a time step At = 0.01 and 201
mesh points are used in computation.

150

 

 

 

 

  

  

 

 

 

 

 
   

 

 

 

 

 

 

 

 

 

 

 

1.0 d g I I . III!—
0.5 - E I I
0.0 . i I I //—I6
0 1 2 3 4 14
(a) r=1.0
1 0 "g 4' """"" | l‘ I 0 ‘ a
0.5 — 5 0.5 -
0-0 I" I I I III—16 0-0 IO
0 1 2 3 4 7 0 5
(c) r=14
1.0 -I‘= .
0.5 — E
I E
0.0 . o
0 5

 

Figure 5.26: Evolution of martensitic fractions at r = 1.0, 1.2, . . . , 2.0. Effective stress
a is normalized by 00. Applied traction is path V, and a time step At = 0.01 and
201 mesh points are used in computation.

151

 

 

 

 

 

 

 

O _ 1 1 . l 1 I 1 1 . 7.

Figure 5.27: Movement of TA (Cf/A interface) and TM (WI/(‘3 interface). Applied
traction is path V, and a time step At = 0.01 and 201 mesh points are used in
computation.

152

 

Chapter 6

Conclusions

The studies in this dissertation deal with the pseudoelastic deformations of shape
memory annular plates subject to traction loads. Using a three—dimensional constitu-
tive model, the responses of austenite plates undergoing forward transformation are
investigated for both the plane stress and the plane strain settings. Further studies
on plates subject to arbitrary traction histories are performed numerically for the
plane stress setting.

The stress-strain relation associated with the martensitic transformation is formulated
based on the J2 deformation theory. The evolution of the martensitic volume fraction
is dictated by the effective stress, and a hysteresis algorithm is used to describe
the history dependent behavior in martensitic transformation. In this model, the
transformation strain is volume preserving, reflecting the predominant character of
martensitic transformation. The hysteresis framework employed can describe the full

transformation behavior as well as the partial transformation behavior.

153

 

First, we studied the plane stress responses of austenite plates subject to forward
transformation. As no reverse transformation is activated and transformation starts
from pure austenite, the martensitic fraction is uniquely determined by the effective
stress. The governing equations are rewritten as a two-variable ordinary differential
equation system on effective stress and an auxiliary variable. The stress fields are
solved by the shooting method. The comparison between the responses of a plate
composed of shape memory material and a plate composed of normal elastic material
shows obvious redistribution of stresses due to phase transformation. The Poisson’s
ratio does not affect the stress responses, and the effective stress is non-increasing
along radial direction. The non-increasing property of the effective stress enables us
to study the FMLS, which is the collection of load pairs just enough to complete
martensitic transformation. Further parametric studies on the FMLS by varying
material constants (Ea, of) and geometry (To/Ti) are performed. Other plots using
the non-increasing property of the effective stress along 7‘ direction are structure
maps, which partition the load in (13,3120) plane into different phase combinations.
The partitioning regions on structure maps for the same material are found to be
dependent on plate geometry. To support a phase distribution in the plate with both
an austenite ring and a martensite ring, the ratio ro/I'I- should be greater than a
critical value.

Then we studied the plane strain responses of austenite plates subject to forward
transformation in a similar fashion. The effective stress and the martensitic faction
are also found to be non-increasing along radial direction. However, the Poisson’s

ratio affects the stress distributions, and numerical examples show that the impact

154

 

of the Poisson’s ratio on the phase transformation in the plate can be significant.
The parametric studies on the F MLs by varying Ea, of, u and ro/rz- are performed.
Again a large Ira/r,- is required for a plate to support an outer austenite ring and an
inner martensite ring.

Finally, we studied the plane stress responses of plates under various continuous
traction histories. The shooting method, which has been successfully used to find the
responses of a global loading austenite plate, is applicable for a subclass of traction
histories that produce either global loading or global unloading responses so that only
one loading type exists on the plate within a time step. However, for traction histories
involving both loading and unloading at the same time, the shooting method can not
attack the problem well. An iterative finite difference program is then developed, and
it solves the plane stress responses of the plates subject to arbitrary traction history
effectively.

By assuming a transformation strain field, the stress can be integrated using the
equilibrium equation and the compatibility equation. The stress-strain relation results
in a nonlinear equation on the transformation strain. The equation is discretized with
a finite difference scheme using an equally-spaced mesh, and is solved numerically
using Newton’s iteration with linear search and backtracking. Numerical results from
such an approach are accurate, and the numerical scheme is robust for solving the
proposed boundary value problem.

Several complicated traction histories are simulated. The results reveal that the re-
sponses of the plates subjected to traction loads are history dependent. The effective

stress along the radial direction need no longer be nonincreasing. Quite often, the

155

 

plate will have loading and unloading concurrently. The evolution of martensitic frac-
tions for different locations on the plate is plotted, which observes the given hystere-
sis model. The appearing, disappearing, and movement of phase interfaces between
austenite and austenite/martensite mixture or martensite and austenite/martensite
mixture for several traction histories are presented.

The structure maps are valuable in finding the phase distribution information espe-
cially when proportional loads are applied to shape memory annular plates. This
information can serve as a general guideline for achieving a desired degree of marten-
sitic transformation in designing shape memory ring devices. The numerical iteration
program can be used to refine the design by obtaining the detailed responses in the
plate under various traction loads.

Some possible further studies include: finding the responses of shape memory plates
subject to temperature changes and traction loads, using multi-variant models to
describing the martensitic transformation behavior, refining the stress-strain relation
to accommodate tension/ compression asymmetry, and investigating the responses of
a plate with displacement boundary conditions.

The transformation behavior will be more complicated if both temperature and
stresses change, especially when the temperature is below martensite starting temper-
ature, the material will exhibit shape memory effect. In these situations, the hysteresis
algorithm employed currently can be extended easily to accommodate phase trans-
formation due to change in temperature and stress, such as the one by Ivshin and
Pence (1994b), where temperature and stress are combined to an equivalent driving

variable. The stress-strain relation needs modification to account for temperature-

156

 

induced transformation strains since J2 deformation theory can not describe non-zero
transformation strains under zero effective stress. One possible solution is to decom-
pose the transformation strain into stress-induced part and temperature-induced part,
with the stress-induced transformation strains following J2 deformation theory and
the temperature-induced transformation strains being isotropic expansion. With such
a change in constitutive model, the analysis of the boundary value problem will be
more difficult.

The two-variant model (Wu and Pence, 1998) can describe the richer transformation
behavior than that of a one-variant model, also this model can treat the temperature-
induced and stress-induced transformation behavior concurrently. This model can be
used for describing the transformation behavior. The three-dimensional constitutive
modeling can then involve superposition of the transformation strains associated with
different martensite variants. Some three-dimensional models (Gillet, 1995, Qidwai
and Lagoudas, 2000) that attempt to reflect tension and compression asymmetry can
also be adopted. The transformation strain components will be described with a
constitutive model involving J3 or 11. Note such a change in constitutive model can
make an analytical study almost impossible. The iterative algorithm develOped can
still be applied as long as the total strain can be decomposed into elastic strain and
transformation strain. However, the evaluation of the Jacobian should be modified
accordingly.

Instead of the traction boundary condition, it is possible to solve the problem for
displacement boundary conditions. For example, if the displacements at r = r,- and

r = r0 are given, then hoop strains at r = 7“,; and r = r0 are known; with the

157

 

constitutive model, two conditions on stresses at these two locations will be obtained
and again we have a two—point boundary value problem. For forward transformation of
an austenite plate, a similar shooting program can be developed to solve the boundary
value problem. For some shape memory ring devices, mixed boundary conditions with
both displacement boundary condition and traction boundary condition may be of
interest. Also shape memory composite ring structures can be analyzed similarly with

appropriate displacement continuity conditions on the interfaces of different materials.

158

Appendix A

Mathematica Code for Numerical
Examples I

A.1 Structure Map for Plane Stress Analysis

(*This is the code for plotting structure map in Figure 3.7*)
(*Initialize Material Constants*)

ri=1.; ro=2.;

sigNORM=75.; sigFTs=75./sigNORM; sigFTf=300./sigNORM;
cE=50000./sigNORM; alpha=0.05; Ea = cE*a1pha;

pi06=Pi/6.0; p103=Pi/3.0; pi02=Pi/2.0; sqrt3=Sqrt[3.0];

(*Definition of g function*)
FTS[s_]:=(s-sigFTs)*Pi/(sigFTf-sigFTs);
FTC = Pi/(sigFTf-sigFTs);

gf[s_]:=Which[s<=sigFTs,O., s>=sigFTf, 1.,True, (1.-Cos[FTS[s]])/2.];
g1f[s_]:=Which[s<=sigFTs,O., s>=sigFTf, 0., True, Sin[FTS[s]]*FTC/2.];

(*Definition of DSigmaDr, DBetaDr*)
Delta[s_,b_]:=Modu1e[{g=gf[s], g1=g1f[s]},
s/2*(1-Ea*(g-s*g1)/(Ea*g+s)*Cos[b+p106]“2)1;
DSigmaDr[s_,b_,r-]:=-s‘2*Cos[b]“2/(r*Delta[s,b]);
DBetaDr[s_,b_,r_]:=Module[{g=gf[sl,g1=g1f[sl}.
s*Cos[b]/(r*De1ta[s,b])*(Sin[b]-Ea*(g-s*g1) /
(2*(Ea*g+s))*Sin[b-p103]) ];

(*Integrate 0DE*)
IntegrateODE[sv_,bv_,rv_]:=Module[{}.
C1ear[r,sF,bF, solnTmp];
solnTmp=NDSolve[{
sF’[r]==DSigmaDr[sF[r],bF[r],r], bF’[r]==DBetaDr[sF[r],bFEr],r],

159

sF[rv]==sv,bF[rv]==bv},{sF[r],bFIrl},{r,ri,ro}];
{{sigma[r-],beta[r_l}}={sF[rI,bFIr]}/.solnTmp;
sigRR[r_]=2./sqrt3*sigma[r]*Sin[beta[r]+pi06];
{sigRRIri],sigRR[ro]} ];

(*Plot Structure Map*)
structureMap=ParametricPlot[
{ IntegrateDDE[sigFTs,bt,ri],IntegrateODE[sigFTs,bt,ro],
IntegrateODE[sigFTf,bt,ro],IntegrateODE[sigFTf,bt,ri] },
{bt,0., 2.*Pi}, Piotaange—>{{-3o,30},{-6,6}} 1;

A2 Shooting Method for Plane Stress Analysis

(*This is the code for plotting Figure 3.10*)
(*Initialize Material constants*)

ri=1.; ro=2.;

sigNORM=75.; sigFTs=75./sigNORM; sigFTf=300./sigNORM;
cE=50000./sigNORM; alpha=0.05; Ea = cE*alpha;

(*Loading Model*)

FTS[s_]:=(s-sigFTs)*Pi/(sigFTf-sigFTs);

FTC = Pi/(sigFTf-sigFTs);

sf[srr_, soo_] := Sqrt[srr‘2-srr*soo+soo“2];
gf[s_]:=Which[s<=sigFTs,O., s>=sigFTf, 1., True, (1.-Cos[FTS[s]])/2.];
g1f[s_]:=Which[s<=sigFTs,0., s>=sigFTf, 0., True, Sin[FTS[s]]*FTC/2.];

(*Governing Equations*)
Af[srr_, soo_] := With[{s
g1 = g1f[sf[srr,soo]]},
4*s“3 +Ea *(3 srr‘2 * g + s g1 *(srr-2 soo)“2 )l;
Bf[srr_, soo_] := With[{s sf[srr,soo], g=gf[sf[srr,soo]],
g1 = g1f[sf[srr,soo]]},
(4*s“3 +Ea *(4 3‘2 * g + (g-s g1) *(2 3‘2 - 38rr soo) ))
*(srr - $00)];
dSerrf[srr_,soo_,r_]
dSooDrf[srr_,soo-,r_]

sffsrr,soo], g=gf[sf[srr,soo]],

-(srr -soo)/r;
Bf[srr, soo] / Af[srr, soo] / r;

(*Elastic Solution*)
sigOElas[r_,pi_,po_] := ((pi*ri‘2-po*ro‘2)-(po-pi)*
ri‘2*ro‘2/r“2)/(ro‘2-ri“2);

(*Integrate 0DE*)

IntegrateODE[srrT_Real, sooT_Rea1, rT_Rea1]:= Modu1e[{},
Clear[sigmaR,sigma0,sigma, solnTmp, er, soF];

160

 

solnTmp=NDSolve[{er’[r]==dSerrf[er[r],soF[r],r],
soF’[r]==dSooDrf[er[r],soF[r],r],erfrT]==srrT,
soFIrT]==sooT},{er[r],soFIrll,{r,ri,ro},
AccuracyGoal->8, PrecisionGoal->8, MaxSteps->1000];

{{sigmaR[r_],sigma0[r-]}}={er[r],soF[r]}/.solnTmp;

sigma[r_Real]=sf[sigmaR[r], sigma0[r]];

sigmaR [ro] ] ;

(*Shooting: Starting From Elastic Solution*)
Shooting[pi_,po_]:=Module[{soTry}, soTry=sig0Elas[ri,pi,po];
FindRoot[IntegrateODE[-pi, sooT, ri]==-po,{sooT,soTry}];
Plot[{sigmaR[r],sigmaDEr],sigmaIr]},{r,ri,ro},PlotRange->A11] ];
Shootingf5.,0.5];

A.3 Structure Map for Plane Strain Analysis

(*This is the code for plotting Figure 4.6*)
(*Initialize Material constants*)

ri=1.; ro=2.;

sigNORM=75.; sigFTs=75./sigNORM; sigFTf=300./sigNORM;
cE=50000./sigNORM; a1pha=0.05; Ea = cE*a1pha; v = 0.3;
pi06=Pi/6.0; pi03=Pi/3.0; piO2=Pi/2.0; sqrt3=Sqrt[3.0];

(*Loading Mode1*)
FTS[s_]:=(s-sigFTs)*Pi/(sigFTf-sigFTs);
FTC = Pi/(sigFTf-sigFTS);

gf[s_]:=Which[s<= sigFTs,0., s>=sigFTf, 1., True, (1.-Cos[FTS[s]])/2.];
g1f[s_]:=Which[s<= sigFTs,0., s>=sigFTf, 0., True, Sin[FTS[s]]*FTC/2.];

(*Governing Equations*)
delta[s_,b_]:=Modu1e[{g=gf[s], g1=g1f[s],sin2bpi06=Sin[2*b+pi06]},
(1-v‘2)*s/2+3/8*Ea“2*g*g1+Ea/4*(1+V)(g+s*g1)+(1-2*v)/4*
Ea*((1+sin2bpio6)*g+(1-sin2bpi06)*s*g1)];
dSigDr[s_,b_,r_]:=Module[{g=gf[s],cosb=Cos[b]},
-cosb“2*((1-v“2)s“2+Ea*(2-v)*s*g+3.*Ea‘2*g‘2/4)/(r*delta[s,b])];
dBetaDr[s_,b_,r_]:=Module[{cosb=Cos[b],sinb=Sin[b],g=gf[s],g1=g1f[s]},
cosh/(4*r)*((4*(1-v‘2)*s+3*Ea*g*(Ea*g1+1)
+(5-4*v)*Ea* s*g1)*sinb+sqrt3*Ea*(1-2v)*(g-s*g1)*cosb)/delta[s,b]l;

(*Integrate 0DE*)

IntegrateODE[8v_,bv_,rv-]:=Module[{},
Clearst,bF,solnTmp];
solnTmp=NDSolve[{sF’[r]==dSigDr[sF[r],bF[r],r],

161

bF’[r]==dBetaDr[sF[r],bF[r],r],sF[rv]==sv,

bF[rv]==bv},{sF[r],bF[r]},{r,ri,ro}];
{{sigma[r_],beta[r_]}}={sF[r],bFIr]}/.solnTmp;
sigmaRBar[r_]=2./sqrt3*sigma[r]*Sin[beta[r]+p106];
sigmaOBar[r_]=2./sqrt3*sigma[r]*Sin[beta[r]-pi06];
xi[r_]=gf[sigma[r]];
sigmaZ[r_]=cE/(1-*v)*(v/cE+alpha*xi[r1/2/sigma1r])

*(sigmaRBar[r]+sigma0Bar[r]);
sigmaR[r-]=sigmaRBar[r]+sigmaZ[r];
sigmaO[r_]=sigmaDBar[r]+sigmaZ[r];
{-sigmaR[ri],-sigmaR[ro]} ];

(*Plot Structure Map*)

structureMap = ParametricPlot[{
IntegrateODE[sigFTs,bt,ri],IntegrateODEIsigFTs,bt,ro],
IntegrateODE[sigFTf,bt,ro],IntegrateODE[sigFTf,bt,ri]},
{bt,0., 2.*Pi}, PlotRange->A11, AspectRatio->Automatic];

A.4 Shooting Method for Plane Strain Analysis

(*This is the code for plotting structure map in Figure 4.9*)
(*Initialize Material Constants*)

ri=1.; ro=2.;

sigNORM=75.; sigFTs=75./sigNORM; sigFTf=300./sigNORM;
cE=50000./sigNORM; a1pha=0.05; Ea = cE*a1pha; v = 0.3;

pi06=Pi/6.0; p103=Pi/3.0; piD2=Pi/2.0; sqrt3=Sqrt[3.0];

(*Loading Mode1*)

FTS[s_]:=(s-sigFTs)*Pi/(sigFTf-sigFTs);

FTC = Pi/(sigFTf-sigFTs);

gf[s_]:=Which[s<=sigFTs,O., s>=sigFTf, 1., True, (1.-Cos[FTS[s]])/2.];
g1f[s_]:=Which[s<=sigFTs,0., s>=sigFTf, 0., True, SinIFTS[s]]*FTC/2.];

(*Governing Equations*)

delta[s_,b_]:=Modu1e[{g=gf[s], g1=g1f[s],sin2bpioS=Sin[2*b+pi06]},
(1-v“2)*s/2+3/8*Ea‘2*g*g1+Ea/4*(1+V)(g+s*g1)+(1-2*v)/4*
Ea*((1+sin2bp106)*g+(1-sin2bp106)*s*g1) ];

dSigDr[s_,b_,r-]:=Modu1e[{g=gf[s],cosb=Cos[b]},
-cosb“2*((1-v‘2)s“2+Ea*(2-v)*s*g+3.*Ea“2*g“2/4)/(r*delta[s,b]) J;

dBetaDr[s_,b_,r_]:=Module[{cosb=Cos[b],sinb=Sin[b],g=gf[s],g1=g1f[s]},
cosb/(4*r)*((4*(1-v“2)*s+3*Ea*g*(Ea*g1+1)+(5-4*v)*Ea* s*g1)*
sinb+sqrt3*Ea*(1-2v)*(g-s*g1)*cosb)/de1ta[s,b] ];

162

(*Elastic Solution & Find {sig, beta} Given {srr, soo} *)
sigRElas[r_,pi_,po_l:=((pi*ri‘2-po*ro‘2)+(po-pi)*ri‘2*ro“2/r“2)/
(ro‘2-ri‘2);
sigOElas[r_,pi_,po_]:=((pi*ri‘2-po*ro“2)-(po-pi)*ri‘2*ro‘2/r“2)/
(ro‘2-ri“2);
sigZElas[r_,pi_,po_]:=v*(sigRElas[r,pi,po]+sig0Elas[r,pi,po]);
sigElas[r_,pi_,po_]:=Modu1e[{srT,soT},srT=sigRElas[r,pi,po];
soT=sig0E1as[r,pi,po]; Sqrt[srT‘2+soT“2-srT*soT]];
getBeta[srr_Rea1, soo_Real, s_Rea1]:=Modu1e[{b, eps=1.0/10“6},
b=ArcSin[sqrt3/2*srr/s]-pi06;
Which[Abs[soo-2/sqrt3*Sin[b-pi06]*s]<eps, betaValue=b,
Abs[soo-2/sqrt3*Sin[2*Pi/3-b-p106]*s]<eps, betaValue=2*Pi/3-b,
True, Abortfll;
betaValue ];

stter[srr_Rea1,soo_Real,sTry_Real]:=Modu1e[{g,szz, s},
s=sTry; g=gf[s];
szz=(soo+srr)*(Ea*g+2 s v)/(2 (Ea*g+s));
Sqrt[((srr-soo)‘2+(soo-szz)‘2+(szz-srr)‘2)/2] ];

getSigmaIsrr_Real,soo_Rea1]:=Module[{szz,soln},
szz=(soo+srr)*v;
sRootInit = Sqrt[((srr-soo)“2+(soo-szz)“2+(szz-srr)‘2)/2];
soln=FindRoot[(stter[srr,soo,sRoot])==sRoot,{sRoot,sRootInit},

AccuracyGoa1->8, PrecisionGoa1->8,MaxIterations->200];

sRoot/.soln ];

getSzz[srr_Real, soo_Real]:=Modu1e[{s,g},
s=getSigma[srr,soo]; =gf[s];
(2*s*v+Ea*g)(srr+soo)/(2*(s+Ea*g)) l;

(*Integrate ODE with serri=srrv and soeri=soov*)
IntegrateODE[srrv_Rea1,soov_Rea1]:=Module[{rv,sv, bv,szz},
rv=ri;
sv=getSigmaIsrrv, soov];
szz=getSzz[srrv, soov];
bv=getBeta[srrv-szz, soov-szz,sv];
Clear[sF,bF,solnTmp];
solnTmp=NDSolve[{sF’[r]==dSigDr[sF[r],bF[r],r],
bF’[r]==dBetaDr[sF[r],bF[r],r],strv]==sv,bF[rv]==bv},

{sF[r],bF[r]}.{r,ri,ro},

AccuracyGoa1->8, PrecisionGoa1->8, MaxSteps->1000];
{{sigma[r_],beta[r_]}}={sF[r],bF[r]}/.solnTmp;
sigmaRBar[r_]=2./sqrt3*sigma[r]*Sinfbeta[r]+p106];
sigmaOBar[r_]=2./sqrt3*sigma[r]*Sin[beta[r]-p106];
xi[r_]=gf[sigma[r]];

163

sigmaZ[r_]=cE/(1-2*v)*(v/cE+alpha*xi[r]/2/sigma[r])*
(sigmaRBar[r]+sigma0Bar[r]);

sigmaR[r-]=sigmaRBarfrl+sigmaZ[r];

sigmaO[r_]=sigma0Barfrl+sigmaZ[r];

sigmaR [ro] ] ;

(*Shooting: Starting From Elastic Solution*)
Shootinngi_,po_]:=Module[{sooInit},
sooInit=sig0Elas[ri,pi,po];
FindRoot[IntegrateODE[-pi,soo]==-po,{soo,sooInit}];
Plot[{sigmaR[r],sigma0[r],sigma[r]},{r,ri,ro},PlotRange->All] ];
Shooting[8.,1.];

A.5 Iterative Finite Difference Method for Plane
Stress Analysis

(*This is the code for plotting Figure 5.26(a)*)

(***********************CONSTANTS INIT*************************)
ri = 1. ; r0 = 2. ; NN = 200;

deltaTime = 0.01 ;

pis = # Sin[2. #] & /0 Range[0, 2 Pi, deltaTime];

pos = # Cos[2. #1/2 & /0 Range[0, 2 Pi, deltaTime];

NStep = Dimensionsfpis][[1]];

ZEROPSGARD = 1.0/10.‘6;

DEBUGFLAG = False;

Protectfri, ro, NN, NStep, ZEROPSGARD];

(*Constants for Global Newton Iteration, see nemerical recipe*)
MAXITS = 200;

TOLF = 0.5/10.‘16; TOLX = 1.0/10.“16;(*Accuracy control*)
TOLMIN = 1.0/10.‘6; STPMX = 100.0; ALF = 1.0/10.“4;
Protect[MAXITS, TOLF, TOLMIN, TOLX, STPMX, ALF];

(***********************DISCRETIZATION*************************)
ri2 ri*ri; r02 = ro*ro;

NN1 = NN + 1; NDIM = 2*NN1;

h = (ro - ri)/NN; h2 = h/2;

rs = Rangefri, ro, h];

r23 = Table[rs[[i]]*rs[[i]l, {1, 1, NN1}];

r33 = Table[rs[[i]]*r23[[i]], Ii, 1, NN1}];

(***********************MATERIAL MODEL*************************)
(*Material Constants*)

164

sigNORM = 100. ;
sigFTs = 300. /sigNORM; sigFTf
sigRTs = 200. /sigNORM; sigRTf
cE = 50. *10“3/sigNORM;

alpha = 0.05; Ea = cE*alpha;
pi03 = Pi/3. ; p106 = Pi/6. ;

400. /sigNORM;
100. /sigNORM;

(*Complete Model*)

FTS[s_] := (s - sigFTs)/(sigFTf - sigFTs)*Pi;
rrc = Pi/(sigFTf - sigFTs);

sf[srr_, soo_] := Sqrt[srr“2 - srr*soo + soo‘2];
RTC = Pi/(sigRTs - sigRTf);

RTS[s_] := (s - sigRTf)/(sigRTs - sigRTf)*Pi;
EPSGARD = 1.0/10.“8;

(*g for both loading and unloading*)
gf[s_, s0_, xi0_] := Which[
s > sigFTf, 1.0, s < sigRTf, 0.0,
(*Loading*)
s >= 80,
If[s <= sigFTs, xiO,
gs = If[sO <= sigFTs, x10,
1. - 2.*(1. - xiO)/Max[(1. + Cos[FTSISOII), EPSGARDII;
1. - (1. - gs)/2. *(1. + Cos[FTSIs]l)].
(*Unloading*)
s < s0,
If[s >= sigRTs, xiO,
gs = If[sO >= sigRTs, xiO,
2. *xiO/Max[(1. - Cos[RTSIsOII), EPSGARDll;
gs/2. *(1. - Cos[RTst]])]
];

(*dg/ds for both loading and unloading*)
g1f[s_, sO_, xi0_] := Which[s > sigFTf, 0.0,
s < sigRTf, 0.0,
(*Loading*)
3 >= 30,
If[s <= sigFTs, 0,
gs = If[sO <= sigFTs, xi0,
1. - 2. *(1. - xiO)/Max[(1. + Cos[FTS[sO]]), EPSGARDl];
(1 gs)*FTC/2. *Sin[FTS[s]]],
(*Unloading*)
s < s0,
If[s >= sigRTs, 0,
gs = If[sO >= sigRTs, xiO,

165

2. *xiO/Max[(1. - Cos[RTS[sOIJ), EPSGARDII;
gs*RTC/2. * 31n[Rrs[s]]]
J;

(***********************ALLOCATE MEMORY*************************)
ZEROVECTOR = TablefO. , {1, 1, NN1}];
ZEROVECTOR2 = Tab1e[0. , {1, 1, NDIM}];
Protect[ZEROVECTOR, ZEROVECTOR2];
InitMatrice[] := Module[{},
(*Integrations computation*)
HS = ZEROVECTOR; Is = ZEROVECTUR;
(*Jacobian related*)
HR = TableIO. , {1, 1, NN1}, {j, 1, NN1}];

H0 = Table[0. , {1, 1, NN1}, {j, 1, NN1}];
IR = Table[0. , {1, 1, NN1}, {j, 1, NN1}];
I0 = Table[0. , {1, 1, NN1}, {j, 1, NN1}];
DSRDER = Table[0. , {1, 1, NN1}, {j, 1, NN1}];
DSRDEO = Table[0. , {1, 1, NN1}, {3, 1, NN1}];
DSODER = TablefO. , {1, 1, NN1}, {j, 1, NN1}];
DSODEO = TablefO. , {1, 1, NN1}, {j, 1, NN1}];

DEDE = Tablefo. , {1, 1, 2* NN1}, {j, 1, 2* NN1}];
(*Working variables*)
SigRl = ZEROVECTUR; Sig01 = ZEROVECTDR; Sigl = ZEROVECTOR;
X11 = ZEROVECTOR; DXiDsl = ZEROVECTOR;
EtrRl = ZEROVECTOR; Etr01 = ZEROVECTOR;
SigO = ZEROVECTOR; XiO = ZEROVECTOR;
(*EtrRO, EtrOO are allocated by Take function*)
1;

(***********************MATRIX; COMPUTE DSDE**********************)
(*HR, HO*)
PrepareDHI] := Module[{i,j},

For[i = 2,1 <= NN1,

HRIIi, 1]] = cE*h2*rs[[1]]:
Ho[[1, 1]] = HR[[1, 1]] + cE*rs[[1]]‘2;
Forfj = 2, j < i,
HRIIi, jll = cE*h*rs[[j]l;
H0[[i, j]] = HR[[1,j]];
j++];
HR[[i, 1]] = cE*h2*rs[[i]];
H0[[i, 1]] = HR[[i, 1]] - cE*rs[[i]]“2;
i++];
];
(*IR, 10*)

166

PrepareDI[] := Module[{i, j, 10, 11},
For[i = 2, i <= NN1,
For[j = 1, j <= i, 10 = Max[j, 2];
For[l = 10, 1 <= i, 11 = l - 1;
IRIIi, jll += h2*(HR[[1, j]]/rs[[1]l“3 +
HR[[11, j]]/rs[[11]]‘3);
I0[[i, jll += h2*(H0[[l, jll/rs[[1]]“3 +
H0[[l1, j]]/rs[[11]l“3);
1++1;
j++];
i++];

1;

(*DSRDER, DSRDEO, DSODER, DSODEO*)
PrepareDSf] := Module[{i, j},
rc1 = rs[[NN1]]‘2/(rs[[NN1]]“2 - rs[[1]]‘2);

For[i = 1, i <= NN1,
rc2 = (rs[[1]]“2/rs[[i]]“2 - 1)*rc1;
rc3 = -(rs[[1]]“2/rs[[i]]‘2 + 1)*rc1;
rc4 = rs[[i]]“2;

For[j = 1, j <= NN1,
DSRDERIIi, 3]] rc2*IR[[NN1, j]] + IRIIi. jll;
DSRDEOffi, jll rc2*IO[[NN1, jll + I0[[i, jll;
DSODERIIi, j]] rc3*IR[[NN1, jl] +
HRIIi. jlllrc4 + IR[[i, jll;
DSODEOIIi, j]] = rc3*IO[[NN1, j]]
HOIIi, jll/rc4 + I0[[i, jll;
j++];
1++];

];

+

PrepareMatrice[] := (InitMatrice[]; PrepareDH[];
PrepareDI[]; PrepareDS[];);

(***********************MATRIX : COMPUTE DEDE**********************)
ComputeHI[] := Module[{hT, iT, i, j, i1},
hT = 0.0; Hs[[1]] = 0.0;
For[i = 2, i <= NN1,
HsIIill = r2s[[1l]*Etr00[[1]] - r23[[i]l*Etr00[[i]];
11 = i - 1;
hT += rsffi]]*(EtrRO[[ill + Etr00[[i]]);
hT += rs[[illl*(EtrRO[[i1]] + EtrOOffilll);
Hs[[1]] += h2*hT;
Hs[[i]] *= cE;
i++];

167

Is[[1]] = 0.0; iT = 0.0;

For[i = 2, i <= NN1,
i1 = i - 1;
1T += Hs[[i]]/r331[i]]; iT += Hsffilll/r3sffilll;
Is[[1]] = 02*1r;

i++];

];

ComputeDEDEII := Module[{i, j, 11, j1},
For[i,= 1, i <= NN1,
For[j = 1, j <= NN1,
ii = i + NN1; j1 = j + NN1;
DEDEIIi, jll = XA[[i]]*DSRDER[[i, jll + XB[[i]]*DSODER[[i, jll;
0000111. 3111 = x111111*0s0000111, 311 + x01111140s0000111, 311;
DEDEIIil, j]] = XBIIi]]*DSRDER[[i, jl] + XCIIi]l*DSODER[[i, j]];
00001111, 3111 = x011111*030000111, jl] + x011111*030000111, jJ];
'++];
1+1];

1;

(***********************MATRIX ; COMPUTE JACOBIAN**********************)
ComputeFVec[] := Module[{i, coeff},
Clear[EtrRO, Etr00];
EtrRO = Take[x, NN1]; Etr00 = TakeIx, -NN1];
ComputeHI[];
c0 = (pi*ri“2 - (po + ISIINNIJI)*ro‘2)/(ro‘2 - ri“2);
c1 = (po - pi + Is[[NN1]])*ri“2*ro“2/(ro“2 - ri“2);
For[i = 1, i <= NN1,
SigleIill co + c1/rs[[1]]“2 + Is[[i]];
Sig01[[i]l CG + (Hsffill - c1)/rs[[i]]“2 + Isffil];
Sig1[[i]] = stSigRiIIill, Sig011[i]]];
Xi1[[i]] = gfISiglffill, SigOIIill, Xi0[[i]]];
0x1Ds1[[1]] = g1f[Sig1[[i]], SigOIIill, XiO[[i]l];
coeff = a1pha*X11[[i]]/2/Sig1[[il];

EtrR1[[i]] = (2 Siglefill - Sig01[[i]l)*coeff;
Etr01[[i]] = (2 Sig01[[i]l - SigleIi]])*coeff;
i++];
fvec = x - JoinfEtrRI, Etr01];

];

ComputeJacobianfl := Module[I},
ComputeFVec[];
XA = TableIalpha/2/Siglffill*(2*Xil[[il] +
(2 SigleIill - Sig01[[i]])‘2/2/Sig1[[il]*(DXiDsl[[ill -
Xi1[[i]]/Sig1[[ill)), {1, 1, NN1}];

168

X8 = TableIalpha/2/Sig1[[i]]*(-Xi1[[i]] +
(2 Siglefill - Sig01[[i]l)*(2 Sig01[[i]] -
SigR1[[i]])/2/Sig1[[i]]*(DXiDsi[[i]] - X11[[1]]/31g1[[i]])).
{1, 1, NN1}];

xc = Tablefalpha/2/Sig1[[i]]*(2*Xi1[[i]] +
(2 Sig01[[i]] - SigR1[[i]])‘2/2/Sig1[[i]]*
(Dx1Ds1[[1]] - X11[[1]]/31g1[[1]])), {1, 1, NN1}];

ComputeDEDEfl:

fjac = IdentityMatrix[2*NN1] - DEDE;

]

(********************GLOBAL NEWTON ITERATION *******************)
LinearSearch[] := Module[{i, j, a, b, 1m, lm2, lmT, lmMin,
rshl, rsh2, tmp, slope, disc},
checkFlag = 0;
pNorm = Normfp];
If[pNorm > stepMax, p = p*(stepMax/pNorm)];
slope = g . p;
If[slope >= 0. ,
Print["Roundoff problem in Linear Search, Abort!"]; Abort[]];
test = 0.0;
For[i = 1, i <= NDIM,
tmp = Abs[p[[i]ll/Max[Abs[xold[[i]]], 1.0];
If[tmp > test, test = tmp];
i++1;
lmMin = TOLX/test;
1m = 1.0;
While[True,
x = xold + lm*p;
ComputeFVeCI]; f = Normffvec]/2.0;
Which[1m < lmMin,
(x = xold; checkFlag = 1;

If[DEBUGFLAG, Print["Convergence on X! Check Root Converge"]];

lmUsed = lm; Break[]),
f <= fold + ALF*1m*slope,
(lmUsed = lm; BreakIJ),
True,
If[lm == 1.0, lmT = -slope/(2.0*(f - fold - slope)),
rhsl = f - fold - lm*slope; rhs2 = f2 - fold - lm2*slope;
a = (rhsl/lm‘2 - rhs2/lm2“2)/(lm - lm2);
b = (-1m2*rhsl/lm‘2 + lm*rhs2/lm2‘2)/(1m - lm2);
If[a == 0.0, lmT = -slope/(2.0*b),
disc = b*b - 3.0*a*slope;
Which[disc < 0.0, lmT = 0.5*1m,
b <= 0.0, lmT = (-b + Sqrt[disc])/(3.0*a),

169

True, lmT = -slope/(b + Sqrtfdiscl);

J;
1;
If[lmT > 0.5*1m, lmT = O.5*lml;
J;
1;
lm2 = lm;
f2 = f;
lm = Max[lmT, 0.1*lm];
1;

J;

GlobalNewton[] := Module[{i, j, test, tmp, den},
ComputeFVec[];
If[Max[Abs[fvec]] < 0.01*TOLF,
Print["Solution found!f=F.F/2=", Norm[fvec]/2.0]; checkFlag = 0;
Goto[return]];
stepMax = STPMX*Max[ Norm[x], NDIM];
f = Norm[fvec]/2.0;
For[iterID = 1, iterID <= MAXITS,
If[DEBUGFLAG, Print["Start of Iterationz", iterIDJ];
ComputeJacobian[];
g = fvec . fjac;
xold = x; fold = f;
p = LinearSolve[fjac, -fvec];
LinearSearch[];
If[MafobsEfveCJJ < TOLF, checkFlag = 0; Goto[return]];
If[DEBUGFLAG,

Print["Result of Linear Search: f=F.F/2=", f, " lm=", lmUsed]];
If[checkFlag != 0, If[DEBUGFLAG, Print["checkFlag not Zero!"]];
den = Max[f, 0.5*NDIM]; test = 0.0;

For[i = 1, i <= NDIM, i++.

tmp Absfgffill*MafobsfxfIilll, 1.0l/den];

If[tmp > test, test = tmp];];

If[test < TOLMIN, checkFlag = 1, checkFlag = O];

Goto[return];

1;

test = 0.0;

For[i = O, i < NDIM, i++,
tmp = Abs[x[[i]] - xoldffilll/Max[Abs[x[[i]ll, 1.0];
If[tmp > test, test = tmp];

J;

If[test < TOLX, Print["\[De1ta]X convergence tested!"l;];

iterID++;

170

Print["MAXITS exceeded in GlobalNewton");
Labe1[return];
If[ f > TOLF, Print["f=", f, " is too large!",
"Decrease step size and retry!"]; Abort[]];
Print["Solution found after ",iterID, " iterations! f=F.F/2=", f];
ReturnEf];
];

(***********************MAIN PROGRAM*********************)
MainProg[] := Module[{i},

(*Initialization*)
PrepareMatrice[];

For[iStep = 1, iStep <= NStep,
If[iStep == 1, x = Table[0.0, {1, 1, NDIM}]];
Unprotectfpi,po];
pi = pis[[iStep]]; po = pos[[iStep]];
Protect[pi, po];
Print["++++++++Step:", iStep, "/", NStep, "+++++++++"];
Print["pi=", pi, " & po=", po];
If[Absfpi] < ZEROPSGARD && Absfpo] < ZEROPSGARD,
(SigRl = ZEROVECTOR; Sig01 = ZERUVECTOR; Sig1 = ZEROVECTOR;
Xi1 = ZEROVECTOR; x = ZEROVECTOR2;
Print["Zero Load! All responses will be zero!"]),
GlobalNewton[]
];
Sigs[iStep] = Sigl; SigRsfiStep] = SigRl; Sig08[iStep] = Sig01;
XisfiStep] = Xi1;
For[i = 1, i <= NN1,
Sig0[[i]] = 31g111111; XiOIIil] = x11[11]];
i++];
iStep++];
1;
MainProg[];

(*Plot Figure 5.26(a)*)
ListPlot[{Sigs[#][[1]]. XiSI#lfflll} & /@ Rangefl. NStepll;

171

Bibliography

Berman, L. B. and White, S. R. (1996). Theoretical modeling of residual and trans-
formational stresses in SMA composites. Smart Structures and Materials, 5:731—743.

Bernardini, D. and Pence, T. J. (2002a). Models for one-variant shape memory mate-
rials based on dissipation functions. International Journal of Non-Linear Mechanics,
37(8):1299—1317.

Bernardini, D. and Pence, T. J. (2002b). Shape-memory materials, modeling. In
Schwartz, M., editor, Encyclopedia of Smart Materials, pages 964-979. John Wiley
& Sons Inc., New York, USA.

Birman, V. (1997). Review of mechanics of shape memory alloy structures. Applied
Mechanics Review, 50:629—645.

Birman, V. (1999). Analysis of an infinite shape memory alloy plate with a circular
hole subjected to biaxial tension. International Journal of Solids and Structures,
36(1):167—178.

Bondaryev, E. N. and Wayman, C. M. (1988). Some stress-strain-temperature rela-
tionships for shape memory alloys. Metallurgical Transactions A, 19Az2407-2413.

Boyd, J. G. and Lagoudas, D. C. ( 1996). A thermodynamical constitutive model for
shape memory materials. Part I. The monolithic shape memory alloy. International
Journal of Plasticity, 12(6):805—842.

Briggs, J. P. and Castaneda, P. P. (2002). Variational estimates for the effective
response of shape memory alloy actuated fiber composites. Journal of Applied Me-
chanics, 69(4):470—480.

Briggs, J. P. and Ostrowski, J. P. (2002). Experimental feedforward and feedback
control of a one-dimensional sma composite. Smart Material and Structure, 11(1):9—
23.

Brinson, L. C. and Lammering, R. (1993). Finite element analysis of the behavior
of shape memory alloys and their applications. International Journal of Solids and
Structures, 30(23):3261—3280.

Budiansky, B. (1958). Extension of Michell’s theorem to problems of plasticity and
creep. Quarterly of Applied Mathematics, 16:307—309.

172

Budiansky, B. (1959). A reassessment of deformation theories of plasticity. Journal
of Applied Mechanics, 26:259—264.

Budiansky, B., Hutchinson, J. W., and Lambropoulos, J. C. (1983). Continuum
theory of dilitant transformation toughening in ceramics. International Journal of
Solids and Structures, 19(4):337—355.

Cash, J. R. (1986). Numerical integration of non-linear two-point boundary value

problems using iterated deferred corrections-I: A survery and comparision of some
one-step formulae. Comp. £5 Maths. with Appls., 12A(10):1029-1048.

Cornelis, I. and Wayman, C. M. (1976). Experiments on hysteresis in a thermoelastic
martensitic transformation. Scripta Metallurgica, 10(4):359—364.

Daele, M. V. and Cash, J. R. (2000). Superconvergent deferred correction methods
for first order systems of nonlinear two-point boundary value problems. SIAM J.
Sci. Comput, 22(5):1697—1716.

Delaey, L., Krishnan, R. V., Tas, H., and Warlimont, H. (1974). Thermoelasticity,
pseudoelasticity and the memory effects associated with martensitic transformations.
Journal of Materials Science, 9:1521—1555.

Duerig, T. W. (2002). The use of superelasticity in modern medicine. MRS Bulletin,
pages 101—104.

Falk, F. (1980). Model free energy, mechanics, and thermodynamics of shape memory
alloys. Acta Metallurgica, 28:1773—1780.

Giannakopoulos, A. E. and Olsson, M. (1993). Axisymmetric deformation of trans-
forming ceramic rings. Mechanics of Materials, 16(3):295—316.

Gillet, Y., Meunier, M. A., Brailovski, V., Trochu, F., Patoor, E., and Berveiller, M.
(1995). Comparison of thermomechanical models for shape memory alloy. Journal
de Physique IV, Colloque, 5(C8):1165—1170.

Huang, M. and Brinson, L. C. (1998). A multivariant model for single crystal shape
memory alloy behavior. J. Mech. Phys. Solids, 46(8):1379—1409.

Huo, Y. and Muller, I. (1993). Nonequilibrium thermodynamics of pseudoelasticity.
Continuum Mechanics and Thermodynamics, 5:163—1993.

Ivshin, Y. and Pence, T. J. (1994a). A constitutive model for hysteretic phase
transition behavior. International Journal of Engineering Science, 32(4):681—704.

Ivshin, Y. and Pence, T. J. (1994b). A thermomechanical model for a one-variant
shape memory material. Journal of Intelligent Material Systems and Structures,

52455—473.

Kamita, T. and Matsuzaki, Y. (1998). One-dimensional pseudoelastic theory of
shape memory alloys. Smart Materials and Structures, 7(4):489—495.

173

 

 

Kuczma, M. S. and Mielke, A. (2000). Influence of hardening and inhomogenuity
on internal loops in pseudoelasticity. Z. Angew. Math. Mech., 80(5):291—306.

Liang, C. and Rogers, C. A. (1990). One-dimensional thermomechanical constitutive
relations for shape memory materials. Journal of Intelligent Material Systems and
Structures, 1:207—228.

Lim, T. J. and McDowell, D. L. (1999). Mechanical behavior of a Ni-Ti shape mem-
ory alloy under axial-torsional proportional and nonprOportional loading. Journal
of Engineering Materials and Technology, 121(1):9—18.

Lim, T. J. and McDowell, D. L. (2002). Cyclic thermomechanical behavior of a
polycrystalline pseudoelastic shape memory alloy. Journal of the Mechanics and
Physics of Solids, 50(3):651—676.

Muller, I. and Xu, H. (1991). On the pseudo-elastic hysteresis. Acta Metall. Mater,
39:263—271.

Nadai, A. (1950). Theory of flow and fracture of solids. McGrawHill, New York.

Novak, V., Sittner, P., Vokoun, D., and Zarubova, N. (1999). On the anisotropy of
martensitic transformations in Cu-based alloys. Materials Science and Engineering
A, 273-275z280—285.

Otsuka, K., Kakeshita, T., and Editors, G. (2002). Science and technology of shape-
memory alloys: New developments. MRS Bulletin, pages 91—100.

Otsuka, K. and Shimizu, K. (1986). Pseudoelasticity and shape memory effects.
International Metals Reviews, 31(3):93—114.

Otsuka, K. and Wayman, C. M. (1998). Shape memory materials. Cambridge
University Press.

Paskal, Y. I. and Monasevich, L. A. (1981). Hysteresis features of the martensitic
transformation of titanium nickelide. Phys. Met. Metall., 52(5):95—99.

Press, W. H., Teukolsky, S. A., Vetterling, W. T., and Flannery, B. P. (2002).
Numerical recipes in 0' ++: The art of scientific computing. Cambridge University
Press.

Qidwai, M. A. and Lagoudas, D. C. (2000). On thermomechanics and transformation
surfaces of polycrystalline NiTi shape memory alloy material. International Journal
of Plasticity, 16(10—11):1309—1343.

Roberts, S. M. and Shipman, J. S. (1972). Two—point boundary value problems:
Shooting methods. American Elsevier Publishing Company, New York.

Saadat, S., Salichs, J., Noori, M., Hon, 2., Davoodi, H., Bar-on, I., Suzuki, Y., and
Masuda, A. (2002). An overview of vibration and seismic applications of NiTi shape
memory alloy. Smart Material and Structure, 11:218—229.

174

 

 

 

Sun, Q. P. and Hwang, K. C. (1993a). Micromechanics modelling for the constitutive
behavior of polycrystalline shape memory alloys—I. Derivation of general relations.
J. Mech. Phys. Solids, 41:1—18.

Sun, Q. P. and Hwang, K. C. (1993b). Micromechanics modelling for the consti-
tutive behavior of polycrystalline shape memory alloys—II. Study of the individual
phenomena. J. Mech. Phys. Solids, 41:19—33.

Tanaka, K. (1990). A phenomenological description on thermomechanical behavior
of shape memory alloys. ASME Journal of Pressure Vessel Technology, 112:158—163.

Tanaka, K., Tobushi, H., and Nishimura, F. ( 1995). Transformation start lines in
TiNi and Fe-based shape memory alloys after incomplete transformations induced
by mechanical and/ or thermal loads. Mechanics of Materials, 19(4):271—280.

Tobushi, H., Tanaka, K., Hori, T., Sawada, T., and Hattori, T. (1993). Pseudoelas-
ticity of TiNi shape memory alloy. JSME Int. J. Series A, 36:314—318.

Trochu, F. and Qian, Y. (1997). Nonlinear finite element simulation of superelastic
shape memory alloy parts. Computers (’9’ Structures, 62(5):799—810.

Wayman, C. M. (1980). Some applications of shape-memory alloys. Journal of
Metals, 32(6):129—137.

Wolfram, S. (2003). The Mathematica Book. Wolfram Media, fifth edition.

Wu, X. and Pence, T. J. (1998). Two variant modeling of shape memory materials:
Unfolding a phase diagram triple point. Journal of Intelligent Material Systems and
Structures, 92335-354.

Yan, W., Wang, C. H., Zhang, X. P., and Mai, Y. (2002). Effect of transformation
volume contraction on the toughness of superelastic shape memory alloys. Smart
Material and Structure, 11:947-955.

175

 

 

 

UNIVERSI LIBRA ES

IIIIIIIIIIIIIIIIIIIIIIII

3 1293 02736 19 5