BIFURCATION PHENOMENON IN BI-LAYERED SYSTEM WITH THIN FILMS AND ITS

APPLICATION IN THE DEVELOPMENT OF STRETCHABLE BENDABLE DEVICES

By

Mayank Sinha

A DISSERTATION

Michigan State University

in partial fulfillment of the requirements

Submitted to

for the degree of

Mechanical Engineering – Doctor of Philosophy

2020

ABSTRACT

BIFURCATION PHENOMENON IN BI-LAYERED SYSTEM WITH THIN FILMS AND ITS

APPLICATION IN THE DEVELOPMENT OF STRETCHABLE BENDABLE DEVICES

By

Mayank Sinha

“Bifurcation of equilibrium” is a term that is used to describe the condition in which a system
admits multiple solutions to the equilibrium problem, in other words has multiple possible deformed
configurations that can be considered equilibrated under the same load, usually referred to as the
critical load. The loss of uniqueness of the equilibrium problem, also results in a loss of stability of
the system, which is because the equilibrated configuration that is unstable becomes energetically
more favorable over the stable one. An unstable equilibrated configuration can be described as a
configuration in which a small change in parameter values, for example a small perturbation in the
load or a small defect in the geometry of the system, will lead to a large change in the behavior of
the system that can often be catastrophic for the system itself, one example of this is the buckling
of a beam.

Conventionally due its disruptive nature, instability is considered undesirable. For this reason,
in the past few decades a lot of studies have focused on modeling and predicting the onset of
instability in many systems and under many loading conditions. Since bifurcation events are
specific for each system, namely are unique after the system is defined, some studies have had the
intention to employ this information to learn more about specific systems. For example, several
studies have focused on understanding how changing the material properties of a system can affect
the onset of instability. The knowledge gained from that can also be employed to solve the inverse
problem – in other words, if we identify the instability onset, we can use this information as a means
to evaluate the mechanical properties of a specific system.

This has proved especially useful for systems in which it is difficult to employ conventional
methods to evaluate the mechanical properties. One example of that, is to use this principle to
estimate mechanical behavior of a bilayered system of silicon nanocrystals (SiNC) deposited on a

layer of polydimethylsiloxane (PDMS). Due to this unique combination of materials, conventional
tests to evaluate mechanical properties, like uniaxial tensile or nano-indentation tests, cannot be
performed. In this thesis, I employed the instability analysis for a bilayered system subjected to
finite bending, to estimate the material behavior of the SiNC layer.

In addition to studying instability in finite bending, I have developed a mathematical framework
to study deformation and instability in a multilayered system undergoing unbending and eversion.
First, I develop the solution of nonlinear elasticity to evaluate the change in geometrical parameters
of the system due to deformation. Then, I solve the instability analysis to evaluate when a bifurcation
occurs for a multi-layered system undergoing unbending or eversion. The solution for the instability
problem is then obtained numerically by employing the compound matrix method.

I have revealed that the stress states in a bilayered structure subject to unbending or eversion
can result in the presence of multiple neutral axes, as previously observed in finite bending as well,
or no neutral axes. Also, I found that the presence of a second layer significantly affects bifurcation
behavior during unbending and eversion. Finally, I show that the presence of the stiff layer can
significantly affect the bifurcation behavior just by moving the stiff layer from outside to inside in
the reference configuration.

TABLE OF CONTENTS

LIST OF TABLES .

. .

LIST OF FIGURES .

. .

.

.

.

.

.

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

vi

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . vii

.

.

.

.

.

.

.

.

Finite bending . .

CHAPTER 1 LITERATURE REVIEW . . . . . . . . . . . . . . . . . . . . . . . . . .

1.2.1
. . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2.2 Unbending and eversion . . . . . . . . . . . . . . . . . . . . . . . . . .

.
1.1 Overview .
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .
1.2 Finite bending and unbending in large deformations . . . . . . . . . . . . . . . . .
.
.
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

1
1
3
3
6
7
. . . . . . . . . . . . . . . 14
1.4.1 Mechanical properties of nanoparticles . . . . . . . . . . . . . . . . . . . . 14
1.4.2
Instability and in situ mechanical properties of thin films . . . . . . . . . . 16
1.5 Objectives of this dissertation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 17

1.3 Bifurcation Theory .
1.4 Silicon Nanocrystals and related mechanical properties

.

.

.

.

.

.

.

2.1.1 Kinematics
2.1.2

CHAPTER 2 NOTATION AND GOVERNING EQUATIONS . . . . . . . . . . . . . . . 19
2.1 Finite Bending of a rectangular block . . . . . . . . . . . . . . . . . . . . . . . . . 19
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 19
Stress Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 21
2.2 Unbending and eversion of a cylinder . . . . . . . . . . . . . . . . . . . . . . . . . 23
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
2.2.1 Kinematics
Stress Description . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
2.2.2
Incremental Bifurcations
. 28
2.3.1 Compound matrix method for a bilayer . . . . . . . . . . . . . . . . . . . . 32

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

2.3

.

.

.

.

.

.

.

.

.

.

.

3.1 Methods .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

CHAPTER 3 USING FINITE BENDING TO EVALUATE PROPERTIES OF SINC . . . 38
. 38
PDMS Fabrication . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1.1
Plasma synthesis of SiNCs . . . . . . . . . . . . . . . . . . . . . . . . . . 38
3.1.2
Silicon Nanocrystal Film Characterization . . . . . . . . . . . . . . . . . . 41
3.1.3
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.1.4 Uniaxial Tensile Test
3.1.5
Finite Bending Experiment . . . . . . . . . . . . . . . . . . . . . . . . . . 42
3.1.6 Modeling bifurcations in elastic layered structures . . . . . . . . . . . . . . 43
3.1.6.1
Finite deformation . . . . . . . . . . . . . . . . . . . . . . . . . 43
3.1.6.2 Bifurcation Analysis . . . . . . . . . . . . . . . . . . . . . . . . 45
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46

3.2 Results and Discussion .

CHAPTER 4 UNBENDING AND EVERSION OF A BILAYER : RESULTS . . . . . . . 50
4.1 Stress distribution in bilayer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.2 Bifurcation behavior during unbending and eversion . . . . . . . . . . . . . . . . . 54
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

4.2.1
4.2.2 Bilayered systems .

Single layer . .

.

iv

4.2.2.1
4.2.2.2

Thinner, stiffer layer at the compressed side (outermost radius) . . 60
. . . . 62
Thinner, stiffer layer at the tensile side (innermost radius)

CHAPTER 5 FUTURE WORK .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

5.1 Evaluating the mechanical properties of thin films deposited on a pre-stretched

. 66

surface .

.

.

. .

. .

. .

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 66

BIBLIOGRAPHY . .

. .

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 68

v

LIST OF TABLES

Table 3.1: Critical angles of bifurcation for different slenderness ratios, expressed in

degrees. Angles are reported for Λ = 2, 3 and 5. . . . . . . . . . . . . . . . . . . 47

vi

LIST OF FIGURES

Figure 1.1: Example of finite bending deformation. Here, the surfaces corresponding to
x = a1 and x = a2 become portions of the curved surfaces of cylinders of
radii r1 and r2 respectively, where r1 > r2. Figure taken from Rivlin (1949) . . .

Figure 1.2: Finite bending of a neo-Hookean three-layer plate shows the presence of
three neutral axes. In this structure, the initial aspect ratio is 1, the ratio of
shear stiffnesses is 20 and ratio of layer thicknesses is 5. Three neutral axes
become visible starting from a bending semi-angle of 56o. To give evidence
of this effect a bending semi-angle of 90o is imposed. Figure taken from
Roccabianca et al. (2010)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.3: A circular cylindrical sector with internal and external radii R1 and R2,
respectively, and sector angle 2Θ0 straightened under plane strain conditions
into a rectangular block of thickness b − a and length 2l. Figure taken from
Destrade et al. (2014a) .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.4: Sketch of bifurcation modes in an elastic incompressible cylinder subject
to uniaxial compression. Here, ‘h’ is the height of the cylinder, and ‘d’
is its diameter. Different bifurcation modes exist for different geometric
configurations of the cylinder. The letters denote different bifurcation modes
and the numbers refer to the circumferential wave number ‘n’. Figure taken
in part from Gei et al. (2004)

. . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.5: To exhibit the phenomenon of plastic buckling: metal sample subject to
compressive load has buckled in the plastic range. Figure taken from Bigoni
(2012)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

.

.

.

.

.

.

.

.

4

5

6

8

9

Figure 1.6: Rock layers subject to longitudinal compressive stresses have bifurcated and

deformed in the plastic regime. Figures taken from Bigoni (2012) . . . . . . . . 10

Figure 1.7: Bifurcation of a two-layer rubber block under finite bending showing the
presence of long wavelength bifurcation modes. Stiffness and thickness ratios
between the layers are (2.687 N/mm2)/(0.095 N/mm2) and (3 mm)/(40 mm).
The stiff layer (86mm × 3 mm × 150 mm, made up of natural rubber, marked
with a white pencil on the sample) is at the compressive side and coats a
neoprene layer (86 mm × 40 mm × 150 mm). Figure taken from Roccabianca
et al. (2010)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 12

. .

. .

.

vii

Figure 1.8: The three scenarios covered by Sigaeva et al. (2018). A cylindrical circular
sector (in the second column), is bent, unbent or everted (in the third column),
and eventually buckles on the face undergoing compression (last column).
Here, κ = αd/αr where αr is the half angle in the reference configuration,
and αd is the half angle in the deformed configuration. Figure taken from
Sigaeva et al. (2018) . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . .

Figure 1.9: A Si particle cluster, deposited on a support wedge, observed by TEM during
multi-stage loading with a W-tip, (a) just prior to contact, (b) holding position
after the first loading stage, (c) after the second loading stage, (d) after the
third loading stage, (e) 0.04 s prior to fracture and (f) fracture during the
fourth loading stage. Figure taken from Lockwood and Inkson (2008) . . . . .

Figure 1.10: The PDMS was stretched uniaxially and then put under a plasma reactor
for SiNC deposition. SEM image showing surface instability formation on
relaxing the prestretched PDMS substrate. Figure taken from Mandal and
Anthony (2016)

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

.

.

.

. 14

. 15

. 17

Figure 2.1: Reference and deformed configurations for the bilayer. The bilayer deforms
from the Cartesian coordinate system in the reference configuration (e1, e2, e3)
to a cylindrical coordinate system in the deformed configuration (er, eθ, ez).
Specifically, a plane at constant x0
1 transforms to a plane at constant r (dashed
line), and a plane at constant x0
2 transforms to a circular arc at constant θ
(dashed-dotted line). Since the out of plane deformation is taken to be zero,
we have x0

.

.

.

.

3 = z.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20

Figure 2.2: Figure showing the unbending deformation of a right circular cylinder to
a cylindrical sector and further into an everted sector. Arrows show the
direction in which the deformation is being modeled . . . . . . . . . . . . . . . 24

Figure 3.1: Steps from plasma deposition to finite bending. SiNCs are deposited through
a slit-shaped orifice onto PDMS, which is then bent to observe formation
of wrinkles on the SiNC surface.
(b)
Schematic diagram of reactor and SiNC synthesis. (c) SiNC/PDMS sample
illuminated with UV light, demonstrating photoluminescence.
(d) Finite
bending applied with a lab-designed apparatus. . . . . . . . . . . . . . . . . . . 39

(a) Photograph of plasma reactor.

Figure 3.2: SEM and FIB cross-sectioning of the SiNC layer. The thickness is estimated
to be (∼4.50 µm). (a) SEM image showing the surface of the SiNC layer on
PDMS, including the FIB-milled area (inset). (b) Higher magnification SEM
image of the FIB-milled area demonstrating the thickness of the SiNC layer.

. . 40

viii

Figure 3.3:

(a) X-ray diffraction confirmed formation of SiNCs. XRD for SiNCs shows
crystal reflections at 2θ values 28.3, 47.3 and 56.1 degrees consistent with
the (111), (220), and (311) planes of silicon.
(b) TEM image confirming
the formation of SiNCs: black circles highlight individual SiNCs. The inset
shows FFT of the correspondig SiNCs of the TEM image.

. . . . . . . . . . . . 41

Figure 3.4: PDMS sample in the finite bending apparatus. (a) Side view, showing sam-
(b) Top view showing

ple (with handles) in the finite bending apparatus.
formation of bifurcations.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43

Figure 3.5: Critical angles of bifurcation evaluated as a function of ratio of material
properties for a PDMS-SiNC bilayer sample. The ratio of the neo-Hookean
0 /µP0 . All curves have been
material parameters is defined as ¯µ0 = µNC
calculated for a thickness ratio H included in the range 0.0013–0.0018. Each
line represent a different value of slenderness ratio, Λ = 2, black solid line; Λ
= 3, black dashed line; Λ = 4, dark gray solid line; Λ = 5, dark gray dashed
line; Λ = 6, light gray solid line. The insert represents a schematic of the
bilayer and defines the slenderness ratio.

. . . . . . . . . . . . . . . . . . . . . 46

Figure 3.6: Plotting points from the experiments on graph obtained by numerical cal-
culations. The ratio of the neo-Hookean material parameters is defined as
0 /µP0 . Shown here are points obtained for Λ = 2, 3 and 5. Insert
¯µ0 = µNC
shows an instability on the surface of a sample with Λ = 2.

. . . . . . . . . . . 47

Figure 4.1: System undergoing unbending and eversion deformation. R(1)

= 2, R(2)

=

e = 4.5, µr = 4. Deformed configurations are for αd = π/2,−π/2. The

4, R(2)
neutral axis for both deformed configurations are as marked . . . . . . . . . . . 51

i

i

i

i

i

i

Figure 4.2: System undergoing unbending and eversion deformation. R(1)

= 2, R(2)

=

4, R(2)
e = 4.5, µr = 20. Deformed configurations are for αd = π/10,−π/2.
The neutral axis for both deformed configurations are as marked. Two neutral
axes observed for deformation with αd = −π/2 . . . . . . . . . . . . . . . . . . 52

Figure 4.3: System undergoing unbending and eversion deformation. R(1)

= 3, R(2)

=

4, R(2)
e = 4.5, µr = 4. Deformed configurations are for αd = π/2,−π/2. The
neutral axis for both deformed configurations are as marked. Two neutral
axes observed for deformation with αd = −π/2 . . . . . . . . . . . . . . . . . . 53

Figure 4.4: Critical deformation angles αd in unbending and eversion, plotted versus
ratio of innermost to outermost radius ratio ρ. Note that, this study was done
for a single-layered system, and that the reference configuration had varying
angles.‘m’ refers to the mode number of bifurcation (Sigaeva et al., 2018)

. . . 54

ix

Figure 4.5: Studying the variation of critical angle of bifurcation θcr with the ratio of
innermost to outermost ratio of the system ρ. Comparing results published
by Sigaeva et al. (2018) for unbending in a single layer, to our own results
obtained by employing compound matrix method on a bilayer. Results for the
single layer were replicated by applying the same material properties to both
layers .

. . . . . . . . . . . . . . . . . . . . . . . . . . . . .

.

.

.

.

.

.

.

.

.

. 56

Figure 4.6: Solving the bifurcation problem to study how instability in different mode
numbers behaves as a function of innnermost to outermost ratio of the system
ρ for a single-layered system undergoing unbending. Solution for mode
numbers from 1 to 12 is presented and the mode number corresponding to
highest critical angle of bifurcation is marked for all ρ to indicate the first
mode which will be observed in a real world deformation. Note the fact that
mode number 1 is the significant mode both at very small, and for higher
values of ρ

.

.

.

.

.

. . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 57

Figure 4.7: Behavior of stretch and angle of deformation for a bilayered system as a
function of innnermost to outermost ratio of the system ρ. Inner layer of the
system has unit thickness, outer layer a thickness of 0.2 and ratio of inner to
outer neo-Hookean coefficients is 10. Similar to Fig.4.6, the mode number
corresponding to highest critical angle of bifurcation is marked for all ρ to
indicate the first mode which will be observed in a real world deformation . .

. 59

Figure 4.8: Variation in bifurcation behavior for bilayered system with outer layer thick-
ness = 0.2, unit thickness for inner layer, and ratio of inner to outer layer
material properties as indicated in the legend. The variation was observed
as a function of the innermost to outermost radius ratio (ρ) for two param-
eters (a) the stretch, and (b) the angle of deformation, at which the system
becomes unstable. The dots indicate the point at which a different mode
number becomes the primary mode of bifurcation . . . . . . . . . . . . . . . . 61

Figure 4.9: Variation in bifurcation behavior for bilayered system with inner layer thick-
ness = 0.2, unit thickness for outer layer, and ratio of inner to outer layer
material properties as indicated in the legend. The variation was observed
as a function of the innermost to outermost radius ratio (ρ) for two param-
eters (a) the stretch, and (b) the angle of deformation, at which the system
becomes unstable. The dots indicate the point at which a different mode
number becomes the primary mode of bifurcation . . . . . . . . . . . . . . .

. 63

Figure 5.1: Schematic showing the process by which samples are made for evaluating
material properties of a SiNC deposited on a pre-stretched substrate. The
substrate, which is initially rectangular, is rolled into a cylinder. The SiNC is
deposited on it using the plasma deposition process outlined in Sec3.1. The
sample is then opened up and bifurcations are observed on the compressed
surface of the deposited layer

. . . . . . . . . . . . . . . . . . . . . . . . . . . 67

x

CHAPTER 1

LITERATURE REVIEW

1.1 Overview

Bodies that undergo large deformations are susceptible to become unstable, and play an impor-
tant role in the failure of structural components due to significant and unpredictable deformations,
often leading to a loss of a system’s load-carrying capacity. Therefore, it follows that prediction of
the onset of instability is an important topic in the mechanics of bodies that undergo and sustain
large deformations. For this reason, much research has focused on this problem in the past few
decades. The earlier works focused on the theoretical analysis of formation of surface instabilities
due to compression of semi-infinite blocks of hyperelastic materials (Biot, 1964). Specifically,
theoretical solutions have been formulated for such instabilities in systems undergoing a variety of
deformation states, such as uniaxial compression, equibiaxial compression, pure shear, finite bend-
ing, and unbending (Triantafyllidis, 1980; Haughton, 1999; Coman and Destrade, 2008; Gent and
Cho, 1999; Destrade et al., 2010, 2014b; Sigaeva et al., 2018; Roccabianca et al., 2010, 2011). Most
of these works focused on mathematically describing the equilibrium bifurcation in homogeneous
media.

Recently, more attention has been focused on the study of instabilities in composite materials,
for example in Roccabianca et al. (2010, 2011) the authors described, both theoretically and
experimentally, the formation of bifurcations in elastic layered structures under finite bending.
These studies suggest that a multi-layered system can behave in a significantly different way than a
homogeneous block of material, and the presence of a stiffness mismatch between the layers plays a
critical role in the identification of the value of the critical load. A very important result by Dryburgh
and Ogden (1999) showed that bifurcations in a multi-laminated system occur at different critical
stretches than for a homogeneous block, and is often dependent on the mechanical and physical
properties of the system. This realization that, in laminated composites, the value of the critical

1

load is directly dependent on the ratio of the mechanical properties of the coupled materials, opened
the door to the development of new techniques that employed the onset of instability for mechanical
characterization. For example, several studies used the formation of surface instabilities to estimate
the material properties of nanocrsyatlline thin films deposited on deformable substrates (Stafford
et al., 2004; Lu et al., 2006; Jiang et al., 2006; Chung et al., 2011). Specifically, after the deposition
of a thin film composed of a single layer or multiple layers of nanocrystals, on the surface of a
soft substrate, like polydimethylsiloxane (PDMS), the composite is deformed until the formation
of surface instabilities. The value of the critical load or the characteristic pattern of the instability
formed, e.g., wavelength, can be used to determine the ratio of mechanical properties between the
thin film and the substrate.

One of the most exciting materials that exists in a nanocrystalline form in silicon. Silicon
is among the most abundant elements in the universe, moreover it is one of the most common
elements on our planet, second only to oxygen by mass. Due to its semiconducting properties, over
the last few decades, silicon has been a major driver of the world’s technological and economical
growth through its applications in technologies like integrated chips and solar cells. However,
there are other interesting properties that silicon exhibits that have not yet been able to be used
in commercially viable technologies. One of these properties is the photoluminescence of silicon
nanocrystals (SiNC) under ultraviolet light. SiNC usually have a diameter smaller than 100 nm,
are inexpensive and easy to to manufacture and exhibit interesting physical and electrical size-
dependent properties. Recent developments have also allowed production of SiNC in the gas phase,
as opposed to solvent-based processes or high temperatures, which facilitates their deposition on
a polymer-based substrate, such as PDMS (Mangolini et al., 2005; Anthony et al., 2011, 2012;
Mandal and Anthony, 2016).

Combining the mechanical properties of a polymer based substrate, with the tunable properties
of SiNC could lead to the introduction of a new class of light-weight and highly bendable / stretchable
technologies, that exhibit the optoelectronic properties of SiNC (Stafford et al., 2006; Park et al.,
2012; Bao and Chen, 2016). However, to be able to design technologies taking advantage of such

2

properties, it is imperative to understand the mechanical characteristics of both the components
individually and together as a system. Previous attempts have been made to understand mechanical
properties of nanocrystals in both the elastic (Greer et al., 2005; Yue et al., 2011; Gaspar et al., 2010)
and plastic (Han et al., 2007; Lockwood and Inkson, 2008; Uchic et al., 2004; Wang et al., 2015,
2017b,a) regimes. However, there still is a significant knowledge gap in the field of mechanical
characterization of thin SiNC film, and this is the focus of this dissertation.

1.2 Finite bending and unbending in large deformations

1.2.1 Finite bending

Finite bending is one of the most commonly found deformation modes for elastic structures, both in
the natural and the technological world. For example, the bending of aortic valves to regulate blood
flow when the heart beats, and the hyponastic response in the leaves of many terrestrial plants. It also
can be observed in a variety of technologies like curved screens in modern day smartphones, steel-
making processes, and in the crafting of many musical instruments. For this reason, the mechanical
description of bodies subjected to finite bending has been studied extensively in the lat 70 years. In
Rivlin (1949), the author described the mechanics of finite bending of an incompressible, isotropic
and nonlinearly elastic rectangular block into a sector of a cylinder. Additionally, the same study
demonstrated the existence of an exact solution, for both bending and straightening (i.e., inverse
kinematics) in the case of a Mooney-Rivlin material description, given by

(cid:16)

3 − 3(cid:17)

2

(cid:16)

WMR =

c1
2

2
1 + λ
λ

2
2 + λ

+

c1
2

2
1λ

2
2 + λ

2
2λ

2
3 + λ

λ

,

(1.1)

where c1 and c2 are material constants and λi (i = 1,2,3) are the principal stretches corresponding
to the deformation. To keep it consistent with linear elasticity, it can be shown that the shear
modulus, µ0 = c1 + c2. Also, it is important to note that when c2 = 0 the Mooney-Rivlin strain
energy function reduces to the neo-Hookean strain energy function

WnH =

c1
2

2
1 + λ

2
2 + λ

µ0
2

2
1 + λ

2
2 + λ

(1.2)

1 − 3(cid:17)

2

2
3λ

3 − 3(cid:1)

2

(cid:0)λ

3 − 3(cid:1) =

2

(cid:0)λ

3

Furthermore, for incompressible materials in plane strain deformation, it can be shown that Mooney-
Rivlin and neo-Hookean descriptions coincide (Bigoni, 2012).

Figure 1.1: Example of finite bending deformation. Here, the surfaces corresponding to x = a1
and x = a2 become portions of the curved surfaces of cylinders of radii r1 and r2 respectively,
where r1 > r2. Figure taken from Rivlin (1949)

Hyperelastic material models, such as neo-Hookean and Mooney-Rivlin, are ideal for describing
the large deformation observed in finite bending in rubber-like materials. A few common examples
of materials in nature that have their material behavior explained by such a model include rubber,
many polymers, and, to a certain extent, soft biological tissue. Since the introduction of neo-
Hookean and Mooney-Rivlin strain energy functions by Rivlin (1948), plenty of other strain energy
functions have been developed for varying applications. Most notable among these include the
Ogden model (for polymers and biological tissue, Ogden (1972)), the Arruda-Boyce model (for
rubber like materials and polymers, Arruda and Boyce (1993)) and the Gent model (for describing
elasticity in rubber, Gent (1996)).

4

Steigmann and Ogden (1997) developed a nonlinear theory for plain strain deformations of an
elastic body which has been coated with a thin elastic film. This theory was applied to analyzing the
equilibrium of a homogeneously deformed half-plane consisting of compressible elastic material
coated along the surface, and numerical results were obtained by employing a compressible neo-
Hookean strain energy function. Later, Dryburgh and Ogden (1999) provided the solution for
plain-strain flexure for the same system as above, but for an incompressible material undergoing a
finite bending deformation. Roccabianca et al. (2010) provided an analytical solution to the plain
strain finite bending of an elastic multilayered thick plate. The solution they provided revealed the
complex stress distributions that are generated inside multilayerd structures for large deformations.
An interesting phenomenon that was observed was the presence of multiple neutral axes in pure
bending (Fig.1.2). Although the presence of multiple neutral axes has been reported before, by Chu
(1998) and Chuang and Lee (2000) (in thermal loading, or due to residual stresses, respectively) this
result is important because it is the first time this has been reported in finite bending deformation.

Figure 1.2: Finite bending of a neo-Hookean three-layer plate shows the presence of three neutral
axes. In this structure, the initial aspect ratio is 1, the ratio of shear stiffnesses is 20 and ratio of
layer thicknesses is 5. Three neutral axes become visible starting from a bending semi-angle of
56o. To give evidence of this effect a bending semi-angle of 90o is imposed. Figure taken from
Roccabianca et al. (2010)

5

1.2.2 Unbending and eversion

In the past few decades there have been a multitude of authors who have developed solutions for
the bending of a rectangular block into a sector of a cylinder, both for the finite deformation and for
the instability analysis, as outlined extensively in section 1.2.1. Few works have been published,
however, focusing on the inverse deformation, namely unbending or straightening of a cylindrical
sector to a rectangular block. Indeed, there have been only a handful of contributions to this subject
(Ericksen, 1954; Aron et al., 1998; Truesdell and Noll, 2004; Aron, 2000, 2005; Destrade et al.,
2014a,b; Sigaeva et al., 2018).

Figure 1.3: A circular cylindrical sector with internal and external radii R1 and R2, respectively, and
sector angle 2Θ0 straightened under plane strain conditions into a rectangular block of thickness
b − a and length 2l. Figure taken from Destrade et al. (2014a)

Ericksen (1954), for the first time, showed that a circular sector can always be straightened to a
rectangular block, by mapping the reference cylindrical polar coordinate system to a rectangular co-
ordinate system (Fig.1.3). Aron et al. (1998), then, discussed the kinematics behind the plain strain
straightening of annular cylindrical sectors composed of an isotropic, compressible and nonlin-
early elastic solid. They employed the compressible neo-Hookean as well as generalized Blatz-Ko
material model to obtain closed-form solutions for the deformation. In (Destrade et al., 2014a,b)
the authors discuss the elastic deformation of a cylindrical sector, composed of an incompressible
isotropic soft solid, that is straightened into a rectangular block. They focused on two different

6

cases of straightening of a sector, involving two sets of boundary conditions, one, of zero moments
at the ends, and the other, of zero forces. The authors showed a significant difference in the stress
state, and also discussed existence, uniqueness, and stability of solution for different strain energy
functions.

Eversion of cylinders, namely the phenomenon in which the inner and outer surfaces of a
cylinder swap roles due to deformation, was studied by Haughton and Orr for both incompressible
(Haughton and Orr, 1995) and compressible (Haughton and Orr, 1997) materials. The neo-Hookean
and Ogden material models were used to evaluate deformations for incompressible materials, and
the Varga and Blatz-Ko material models were employed for the case of compressible materials.
Finally, Sigaeva et al. (2018) developed a solution which could be used in the case of eversion, but
also in the general case of unbending from a cylindrical sector to another cylindrical sector with a
bigger with the same orientation (i.e., decrease in the bending angle).

1.3 Bifurcation Theory

Elastic material under compression loads might develop instabilities, which mathematically are
identified as a bifurcation of the solution of the equilibrium equations. The earliest study on this
phenomenons has been performed by Euler in the 18th century, in which he studied the buckling
of cylindrical columns under a compression axial load. Employing linear elasticity theory, the
author derived the was named after him, the so called Euler’s formula, that gives the maximum
compression load in the axial direction that a slender ideal column can carry without buckling. Since
then, a lot of work has been done on various different systems under compression, such as plates,
torsional buckling of a beam, and compression of cylindrical blocks. Gei et al. (2004) performed
bifurcation analysis of silicon nitride cylinders undergoing uniaxial compression, to study the
modes of bifurcation on a 3-dimensional system. The authors performed experiments on cylinders
with varying aspect ratios (i.e., the height divided by the diameter of the cylinder), and performed
theoretical analysis to find bifurcation modes corresponding to the amount of compression load
applied.

7

Figure 1.4: Sketch of bifurcation modes in an elastic incompressible cylinder subject to uniaxial
compression. Here, ‘h’ is the height of the cylinder, and ‘d’ is its diameter. Different bifurcation
modes exist for different geometric configurations of the cylinder. The letters denote different
bifurcation modes and the numbers refer to the circumferential wave number ‘n’. Figure taken in
part from Gei et al. (2004)

Depending on the material, owing to the large deformations associated with buckling, bifur-
cations may involve irreversible deformations, namely plastic deformations. Fig.1.5 shows a steel
specimen which has been permanently deformed after loading in compression. In nature also we
find plenty of examples of cases which exhibit plastic deformation associated with buckling. For
example, Fig.1.6 shows the bifurcation of rock layers which have been subjected to longitudinal
compressive stresses and have deformed in the plastic regime Bigoni (2012). Failure due to buck-
ling can also be found in man-made systems, for example when pavements made of concrete are
heated by the sun and expand, forcing adjacent pieces to push against one another, they can show
cracking and lifting up.

Buckling of columns represent what we call ”global modes” of bifurcation. There are other

8

Figure 1.5: To exhibit the phenomenon of plastic buckling: metal sample subject to compressive
load has buckled in the plastic range. Figure taken from Bigoni (2012)

type of instability that have a localized effect, such as surface instabilities. One of the first re-
searchers to investigate the phenomenon of surface instability was Biot (1963; 1964). His work
was expanded in Triantafyllidis (1980), where the authors employed incrementally linear constitu-
tive models for studying the bifurcation problem of an incompressible plate under pure bending.
They also performed a FEM calculation to verify the findings of the theoretical formulation. Only
later, for example due to the work done by Steigmann and Ogden (1997), Dryburgh and Ogden
(1999), Haughton (1999), Gent and Cho (1999) and Coman and Destrade (2008), did solutions
for instability formation in nonlinear, isotropic elastic plates become well known. Gent and Cho
(1999), specifically, performed experiments to validate the results from theoretical predictions in
addition to outlining the theoretical calculations for three different cases, (a) simple compression
in one direction, (b) compression in one direction, with perpendicular direction constrained (pure
shear), and (c) equibiaxial compression. In this study, the authors evaluated that surface instability
is expected to occur, for all cases, at a moderate compression strain ranging from 33% to 55%.
Specifically, for unidirectional compression the critical compression strain for surface instability

9

(a) Gently folded rock layers at Constitution Hill,
Aberystwyth, S. Wales, UK. The folding has been the
consequence of buckling on initially straight layers

(b) Severe folding of metamorphic rock layers initiated
as buckling due to compression stress in Trearddur
Bay, Holyhead, N. Wales, UK. The one-pound coin is
present to show scale

Figure 1.6: Rock layers subject to longitudinal compressive stresses have bifurcated and deformed
in the plastic regime. Figures taken from Bigoni (2012)

10

was evaluated to be about 46%. It is important to notice that these solutions where estimated under
the assumption of a infinite half-space under compression. However, experiments revealed that
sharp creases started occurring for a star significantly lower than that, about 35%. An explanation
for this phenomenon was reported in Coman and Destrade (2008), where the authors developed a
method to theoretically estimate the critical strain of bifurcation for a finite size block (as opposed
to an infinite half-space assumption). Defining η to be the ratio of thickness to the length of the
block, they observed that an Euler-type buckling instability is observed for all values of η (i.e.,
0 < η < ∞). Furthermore, the buckling-like instability was shown to degenerate into a surface
instability in the limit η → ∞. Finally, another work published by Destrade et al. (2010) focusing
on the modeling of rubber as a bimodular material, with it being stiffer in compression than in
tension. In this work, it is shown that flexure bifurcations happen at a lower compressive strain
compared to what was predicted by previous papers, which better approximated what observed
experimentally in Gent and Cho (1999).

In 1997, Steigmann and Ogden developed a mathematical formulation to describe the kinemat-
ics of elastic solids coated by thin elastic films, as described in Sec. 1.2.1. In addition, the authors
developed incremental equations that described a small departure from an equilibrium configura-
tion, to investigate the stability of the deformed configuration, and the possibility of bifurcation. On
comparing the influence of coating on the bifurcation behavior of the half-plane with an uncoated
half-plane, it was found that the presence of a coating has a significant effect on critical stretch of
the half plane, depending on the stiffness of the coating. Specifically, for a very stiff coating, the
critical stretch is very close to one (i.e. almost undeformed state), while for a very soft coating
it approaches the critical stretch of the substrate. Dryburgh and Ogden (1999) later solved the
incremental problem for a thin coating on an elastic surface undergoing a non-homogeneous plane
strain deformation, namely finite bending. The results showed that bifurcations occurred in coated
structures at a critical strain lower when compared to uncoated structures. This effect increases
when the shear modulus of the coating increases with respect to that of the substrate.

Roccabianca et al. (2010) developed a mathematical formulation to evaluate the formation

11

Figure 1.7: Bifurcation of a two-layer rubber block under finite bending showing the presence of
long wavelength bifurcation modes. Stiffness and thickness ratios between the layers are (2.687
N/mm2)/(0.095 N/mm2) and (3 mm)/(40 mm). The stiff layer (86mm × 3 mm × 150 mm, made up
of natural rubber, marked with a white pencil on the sample) is at the compressive side and coats a
neoprene layer (86 mm × 40 mm × 150 mm). Figure taken from Roccabianca et al. (2010)

of bifurcations on elastic multilayered structures subjected to finite bending. In their theoretical
analysis, the authors considered the deformation of a stress-free multi-laminated rectangular block,
described in the Cartesian coordinate, into a sector of a cylindrical tube, described by semi-angle ¯θ;
the material of all layers was described by a Mooney-Rivlin strain energy function. They used the
determinant method to evaluate the critical angle of bifurcation, ¯θcr and concluded that, for a thick
substrate, even a thin coating can have a significant effect on the critical angle of bifurcation in pure
bending. In addition to providing a theoretical and numerical basis for their work, experiments
were also designed and performed to substantiate the theoretical findings. These experiments quite
effectively demonstrated that the theory can be very successfully used as a tool for predicting
the onset of bifurcation in multi-layered structure. In a follow up publication Roccabianca et al.
(2011) the compound matrix method was employed (Backus and Gilbert, 1967; Ng and Reid,
1979a,b, 1985) to understand the formation of long-wavelength modes. It was concluded that it is
possible to find and design systems that exhibit long-wavelength bifurcation modes (longer than
the short-wavelength modes present for surface instabilities).

12

There have been few publications focusing on the instability analysis of sectors of cylindrical
tubes straightened into rectangular block (Destrade et al., 2014a,b), or everted into sectors of
cylindrical tube where the inner and outer faces swap roles (literally turned inside out) (Sigaeva
et al., 2018). Specifically, Haughton and Orr (1995, 1997) performed the bifurcation analysis for the
eversion problem for both incompressible and compressible materials. An important result was that
thin walled tubes can be everted into a cylindrical shape, in both compressible and incompressible
materials, while sectors of tubes with a thickness to radius ratio above a certain critical value, will
become unstable upon eversion. Destrade et al. (2014a,b) investigated the elastic deformation of
a cylindrical sector straightened into a rectangular block. In these works, the authors considered
an incompressible and isotropic elastic material. Within the span of two papers they investigated
the incremental boundary value problem, the system of forces required to sustain said deformation,
and the influence of geometric parameters and constitutive models on the appearance of wrinkles
on the compressed surface of the block. The authors employed two different numerical methods
to solve the system numerically, the compound matrix method and the impedance matrix method.
The impedance matrix method resulted being more robust and precise than the compound matrix
method. Building on that, in 2018, Sigaeva et al. generalized the straightening problem to bending,
unbending, and eversion. They reveal that when the body is modeled as a homogeneous, isotropic,
incompressible, and hyperelastic material, then the solution to the problem always is unique and
exists, for a strain energy function satisfying the strong ellipticity condition. They used the Mooney-
Rivlin strain energy function, and performed theoretical and numerical calculations to predict the
number and wavelength of wrinkles forming on the compressed surface. They also outlined explicit
asymptotic solutions for instability in thin sectors. Later they also performed finite element analysis
to simulate the formation of creases on the surface.

13

Figure 1.8: The three scenarios covered by Sigaeva et al. (2018). A cylindrical circular sector (in
the second column), is bent, unbent or everted (in the third column), and eventually buckles on the
face undergoing compression (last column). Here, κ = αd/αr where αr is the half angle in the
reference configuration, and αd is the half angle in the deformed configuration. Figure taken from
Sigaeva et al. (2018)

1.4 Silicon Nanocrystals and related mechanical properties

1.4.1 Mechanical properties of nanoparticles

Nanoparticles are objects with at least one dimension smaller than 100 nm. Due to their size and
structure, nanoparticles can exhibit physical properties which can be very different from their bulk
counterparts. For example, silicon in its nanocrystalline form exhibits photoluminescence which is
not observed in bulk silicon. As such, nanoparticles have a variety of different uses ranging from
light weight composites for use in vehicles and sports equipment, to cosmetic and skin care, to
additives to clean up water and eliminate harmful bacteria.

In the past few decades, many studies have focused on evaluating the mechanical properties of
nanocrystalline materials / nanoparticles. Most of these studies were concerned with estimating the

14

Figure 1.9: A Si particle cluster, deposited on a support wedge, observed by TEM during multi-
stage loading with a W-tip, (a) just prior to contact, (b) holding position after the first loading stage,
(c) after the second loading stage, (d) after the third loading stage, (e) 0.04 s prior to fracture and
(f) fracture during the fourth loading stage. Figure taken from Lockwood and Inkson (2008)

deformation of metallic samples, at the nanoscale, under tension or compression. The deformation
field at the nanoscale has been measured using a scanning electron microscope (SEM, Greer et al.
(2005); Uchic et al. (2004)) or a tunneling electron microscope (TEM, Wang et al. (2015, 2017b)).
Although these experiments provide reasonably accurate result, one major disadvantage for all the
methods mentioned above is that these tests are inherently irreversible and destructive, due to the
kind of plastic deformations being applied to the samples.

For silicon in the nanoscale, Lockwood and Inkson (2008) performed real-time imaging of
nanoindentation on 50 nm Si nanoparticle clusters, placed on a hard support wedge, with a W-tip
indenter. At higher strains the sample fractured at a weak interface between two nanoparticles
within the cluster (Fig.1.9). Fedorchenko et al. (2009) calculated the Young’s modulus for Si
nanofilms deposited on a hard surface. The films were of varying thickness (10-80 nm) and it
was observed that the calculated modulus approached the value for bulk silicon (150 GPa) as the

15

thickness increased. However, none of the aforementioned papers deal with the presence of a soft
substrate where in situ calculations could reveal important contributions to the interaction of these
nanoparticles with a soft substrate.

1.4.2 Instability and in situ mechanical properties of thin films

To evaluate the mechanical property of thin films, methods that employs the onset of instability
formation on the surface of a nanocrsytalline thin film have been developed in the last 20 years.
Among these, inspired in part by the work of Bowden et al. (1998), the most promising method
is the one developed by Stafford et al. (2004), called “Strain-induced elastic buckling instability
for mechanical measurements” (or SIEBIMM, in short). Briefly, a thin film is transferred on to
a thick substrate like PDMS, then the sample is stretched uniaxially and an instability is induced
due to the compression in the lateral direction, because of the Poisson effect. Since then, a lot
of other groups have used SIEBIMM (Nolte et al., 2005; Jiang et al., 2006; Lu et al., 2006) to
evaluate elastic properties of thin films. For example, Stafford et al. (2006) used this method to
evaluate the mechanical properties of polystyrene and polymethylmethacrylate (PMMA) films of
thicknesses ranging from 5 nm to 200 nm and Vendamme et al. (2007) did the same for ultrathin
films of an elastomeric polyacrylate network inter-penetrated by a silica network. Likewise, Torres
et al. (2009) used SIEBIMM to evaluate the elastic modulus of amorphous methacrylate thin films
and their relation with the glass transition temperature for polymers, whereas Tahk et al. (2009)
used it for the evaluating mechanical moduli of organic electronic materials.

These methods to calculate the mechanical properties of nanocrystalline thin films, however,
employ consistently the use of homogeneous deformation (uniaxial tension or compression). Fur-
thermore, the requirement for specialty equipment, like microscopes, inhibit the widespread use of
this method. Finally, this method can only describe the linear regime, i.e., Young’s modulus, and
gives no information on the nonlinear behavior of the thin film.

As mentioned previously, a lot of work has been done in the development of the theory of
instabilities in non-homogeneous deformations like finite bending (Roccabianca et al., 2010, 2011).

16

Figure 1.10: The PDMS was stretched uniaxially and then put under a plasma reactor for SiNC
deposition. SEM image showing surface instability formation on relaxing the prestretched PDMS
substrate. Figure taken from Mandal and Anthony (2016)

A recent study by Mandal and Anthony (2016) has demonstrated the formation of instabilities on a
thin layer of plasma deposited SiNC on a prestretched elastomeric substrate (PDMS). Combining the
onset of instability with the predictions from theoretical modeling of bifurcation in finite bending,
it is possible to evaluate the mechanical properties of SiNC thin film. The methods we will be
proposing are relatively general and could be used for the in situ measurement of the mechanical
properties for a variety of bilayered systems

1.5 Objectives of this dissertation

Taking into account the literature survey performed and the need for an appropriate method for
evaluating the mechanical properties of silicon nanocrystals, the objectives of this dissertation are
as follows

1. Measuring the mechanical properties of thin film SiNC using the formation of instability in

finite bending;

2. Developing the mathematical model for unbending of a bilayered system from a cylinder;

3. Assessing the stability of the system in unbending, and study the appearance of instabilities

due to the deformation;

17

4. Evaluating the mechanical properties of thin film SiNC in an unbending deformation to

validate the results obtained from the finite bending method.

18

CHAPTER 2

NOTATION AND GOVERNING EQUATIONS

2.1 Finite Bending of a rectangular block

Let x

0 denote the position of a material point in the stress-free reference configuration of an
elastic body B0. A deformation is applied to the system such that positions of material points,x
in the current configuration B, are a function of the material point in the reference configuration
given by x = ψ(x0). We define the deformation gradient as F = Grad(ψ). Further we can define
the left and right Cauchy-Green tensors as, B = FFT and C = FTF, respectively. We can define
the principal stretches, λi (i = r, θ, z) for a given deformation as the square root of the eigenvalues
of B and C.

For an isotropic incompressible and elastic material, the constitutive equations can be written

as a relation between the Cauchy stress T and principal stretches λi as follows

Ti = −π + λi

∂W
∂λi

(2.1)

where π is the Lagrangian multiplier, which represents the hydrostatic pressure associated with
incompressibility, and W is the strain energy function associated with the deformation. Also, in this
section we will only consider incompressible materials, therefore det(F) = 1. Since the determinant
of a tensor is equal to the product of its eigenvalues, we get

λr λθ λz = 1

(2.2)

2.1.1 Kinematics

In this section, first a model for the finite bending of multilayered elastic plates is given. Then
the problem of instability is studied by solving the incremental equilibrium equation. The finite
bending deformation is prescribed as follows (see Fig.2.1)

0(s)
• The plane at constant x
1

transforms to a circular arc at constant r(s)

19

Figure 2.1: Reference and deformed configurations for the bilayer. The bilayer deforms from the
Cartesian coordinate system in the reference configuration (e1, e2, e3) to a cylindrical coordinate
system in the deformed configuration (er, eθ, ez). Specifically, a plane at constant x0
1 transforms to
a plane at constant r (dashed line), and a plane at constant x0
2 transforms to a circular arc at constant
θ (dashed-dotted line). Since the out of plane deformation is taken to be zero, we have x0

3 = z.

0(s)
• The plane at constant x
2
• x

0(s)
3

= z(s), since the out of plane deformation is zero

transforms to a circular arc at constant θ(s)

where the superscript s (= 1,2), refer to the two different layers in the system. If h(s) is the current
thickness of the cylindrical section, then the incompressibility constraint imposes that

r(s)

i

=

l0h(s)
2 ¯θh(s) − h(s)
0

2 .

(2.3)

The layers are assumed to be perfectly bonded to each other such that r(s)

+ h(s). Due
to this all, radial coordinates r(s) can be referred to the same origin O of a cylindrical coordinate
system Orθz, which is common to all deformed layers. Therefore, the index s is omitted from the

= r(s−1)

i

i

20

description of all coordinates, and the deformed configuration will be described in terms of the
global system Orθz.

Imposition of the incompressibility constraint and of the boundary conditions at x0
1 = −h(s)

0 /2, provides all the stretches in the multilayer as

and x(s)

2 = ±l0/2

λθ = λ =

2 ¯θr
l0

,

λr =

=

1
λ

l0
2 ¯θr

,

λz = 1,

and the current thickness of the s-th layer, h(s), as a function of h(s−1), is given as
l0h(s)
0
¯θ

(cid:118)(cid:116)(cid:18) l0h(s−1)
2 ¯θh(s−1) + h(s−1)/2

2 ¯θh(s−1) − h(s−1)

h(s) = − l0h(s−1)

(cid:19)2

2

+

+

0

0

.

(2.4)

(2.5)

Therefore, once the current thickness of the first layer, h(1), is known, the current thicknesses

of all other layers can be found out.

2.1.2 Stress Description

The Cauchy stress tensor in the layer s can be represented in polar coordinates as

T(s) = T(s)

r er ⊗ er + T(s)

θ eθ ⊗ eθ + T(s)
and can be calculated from the constitutive equation Eq.(2.1).
In this study, we adopted the
Mooney-Rivlin strain energy function, which can be written, for each layer s, in terms of the
principal stretches as

z ez ⊗ ez

(2.6)

(cid:16)

W(s) =

c(s)
1
2

2
r + λ
λ

2
θ

+ λ

z − 3(cid:17)

2

(cid:16)

c(s)
2
2

+

r − 3(cid:17)

2

2
λ
r λ

2
θ

+ λ

2
θ λ

2
z + λ

2
z λ

,

(2.7)

1 and c(s)

where c(s)
2 are material constants that can be determined by experiments. For consistency
with linear elasticity for small strains, we can define the linearized shear modulus µ0 for a given
layer s as

For plane strain, the definition of the strain energy, Eq.(2.7), simplifies to

0 = c(s)
(s)
1 + c(s)
µ
2 .
(cid:18)

2 +

λ

(s)
µ
0
2

1
2 − 2

λ

.

(cid:19)

ˆW(s) =

21

(2.8)

(2.9)

where, ˆW(s) = W(s)(1/λ, λ,1), ()(cid:48) denotes differentiation with respect to stretch λ. The integration
of the equilibrium equation within each layer yields the following components of Cauchy stress in
terms of principal stretches (simplified from Eq. (2.4))

(cid:18)

(cid:19)

T(s)
r = ˆW(s) + γ
T(s)

λ ˆW(s)(cid:17)(cid:48)

(cid:16)

=

θ

(s) =

+ γ

(s)
µ
0
2
(s) =

(cid:18)

2 +
λ
(s)
µ
0
2

λ

1
+ γ
2
2 − 1
2

3λ

λ

(s)

(cid:19)

(s)

,

+ γ

,

(2.10)

where, and γ(s) is an integration constant.

The constants of integration γ(s)(s = NC, P) and height h(P) can be calculated by imposing the
boundary conditions on the external surfaces and the interfacial continuity conditions. Absence of
forces on the external surfaces of the multilayer yields the following two boundary conditions

r

T(P)
T(NC)

(r(P)
) = 0,
(r(NC)

i

e

r

) = 0,

(2.11)

e

= r(P)

where r(NC)
+ h(P), because of the interface condition ensuring no relative displacement
between the layers. The interface continuity conditions also impose continuity of tractions at
between layers. This provides an equation which can be written as

i

T(P)

r

(r(P)

e

) = T(NC)

r

(r(NC)

i

).

(2.12)

The constant γ(NC) can be calculated from the boundary condition shown in Eq.(2.11)2, and

using the radial component of Cauchy stress as defined in Eq. (2.10)1), which results in

(NC)
(NC) = − µ
0
2

γ

αr(NC)

e

+

(cid:17)2

(cid:35)

(cid:17)2

(cid:16)

1
αr(NC)

e

(2.13)

where, α = 2 ¯θ/l0. Using above equation, Eq. (2.13), and the interface continuity condition, Eq.
(2.12), one can compute the remaining integration constant as

(P) =

γ

(NC)
(P)
− µ
µ
0
0
2

(NV)

+ γ

(2.14)

(cid:17)2

+

(cid:35)

(cid:17)2

(cid:16)

1
αr(P)

e

αr(P)

e

22

(cid:34)(cid:16)

(cid:34)(cid:16)

= r(P)

where r(NC)
solved to determine the value of γ(P).

e

i

. This equation, along with the value of γ(NC) obtained previously, can now be

On simplifying the previous equations (Eqs. (2.13),(2.14)) we get

(cid:34)(cid:0)αr(1)

i

(cid:1)2

(1)
µ
0
2

(cid:35)

(cid:0)αr(1)

1

(cid:1)2

i

+

(1)
(2)
0 − µ
µ
0
2

+

(cid:34)(cid:0)αr(1)

e

(cid:1)2

(cid:0)αr(1)

1

e

(cid:1)2

+

(cid:35)

(2) = 0,

+ γ

(2.15)

i

e

,r(1)

where r(1)
and γ(2) are functions of h(1) from Eqs. (2.3), (2.5), and (2.14). Therefore the
solution to a bilayer subject to finite bending can be determined after the thickness of the first layer
h(1) is obtained by numerically solving the previous equation.
2.2 Unbending and eversion of a cylinder

2.2.1 Kinematics

Consider a bilayered right, circular, hollow cylinder, the undeformed configuration of which is
defined as

R(s)
i ≤ R ≤ R(s)

e , −π ≤ Θ ≤ π,

and 0 ≤ Z ≤ L,

i

e

and R(s)

where (R,Θ,Z) are the cylindrical coordinates, with orthonormal basis (ER, EΘ, EZ), the superscript
s (= 1,2) denotes the layers of the system, L is the axial length, and R(s)
are the inner
and outer radii of a given layer, respectively, for both layers (s = 1,2). As can be seen in Fig.
2.2, the layer labelled (1) is always the inner laye rin the reference configuration. The reference
configuration is taken to be a hollow cylinder, and it extends from an angle of −π to π. The
deformed configuration is achieved by applying end forces and moments on the open radial surface
of the system. Like the reference configuration, the deformed configuration is also in a cylindrical
coordinate system (r,θ,z), with orthonormal basis (er, eθ, ez), and the deformation is described by
the deformed angle, 2αd (see Fig. 2.2). The radii in the deformed configuration are given by r(s)
and r(s)

e , and current axial length is given by (cid:96).
For our case we exclude the possibility of scenarios in which the deformed angle αd is exactly
zero, which can be evaluated only using Cartesian coordinates. Therefore, the deformed angles

i

23

must be included in αd ∈ [−π, π]\{0}.
The deformation can be modeled as

r = r(R),

θ = κΘ ,and

z = λZ,

(2.16)

where λz = (cid:96)/L is the axial stretch, and κ is a measure of the change in the angles. It defined as
the ratio of the angle in the deformed configuration to the angle in the reference configuration. It
follows that

κ =

αd
π

∈ [−1,1]\{0}.

(2.17)

Hence, 0 < κ < 1 corresponds to unbending deformation and κ < 0 corresponds to deformation
beyond straightening, namely the eversion of the cylinder.

Figure 2.2: Figure showing the unbending deformation of a right circular cylinder to a cylindrical
sector and further into an everted sector. Arrows show the direction in which the deformation is
being modeled

Defining the radii in the current configuration as r(s)

e ), in the case
of an unbending deformation (i.e., κ > 0), the deformed cylindrical sector occupies the following
space

= f(R(s)

) and r(s)

e = f(R(s)

i

i

r(s)
i ≤ r ≤ r(s)

e , −αd ≤ θ ≤ αd,

and 0 ≤ z ≤ (cid:96).

24

That is because the faces identified as inner and outer in the reference configuration, are still
identified as inner and outer in the deformed configuration, respectively. On the other hand, when
κ < 0, which corresponds to the case of eversion, the inner and outer faces swap roles in the
deformed configuration with respect to the reference configuration, namely the cylinder is turned
inside out. We can therefore write

r(s)
e ≤ r ≤ r(s)

i

, αd ≤ θ ≤ −αd

and 0 ≤ z ≤ (cid:96)

For this unbending deformation we define the deformation gradient, F as follows

F =

κr

R eθ ⊗ EΘ + λz ez ⊗ EZ

(2.18)

Since we are only focusing on incompressible solids, we can apply the incompressibility constraint
as det(F) = 1, which lead to the following relation

dr

dR er ⊗ ER +
(cid:16)

(cid:16)

(cid:16)

r(s)(cid:17)2
(cid:17)2

r(s)

e

=

=

(cid:17)2

(cid:17)2

R(s)

i

R(s)

i

R(s)(cid:17)2 −(cid:16)
(cid:17)2 −(cid:16)

R(s)

e

(cid:16)

κλz

(cid:16)

(cid:16)

+

r(s)

i

+

r(s)

i

(cid:17)2

(cid:17)2

.

.

κλz
When this relationship is specified for r(s) = r(s)

e , it allows us to calculate r(s)

e = f(R(s)

e ), as

(2.19)

(2.20)

The principal stretches in any deformation can be calculated as the square root of the eigenvalues of
the Left Cauchy-Green tensor B = FFT. Therefore, from Eq. (2.18) we get the following principal
stretches

λr =

R

|κ| λzr ,

λθ =

|κ| r
R ,

and λz.

(2.21)

One thing to note here is that we have employed the use of the absolute value of κ to ensure that
the values of stretch are always positive. This is especially important for cases in which the system
undergoes eversion because for those cases αd < 0, and therefore κ < 0.

25

2.2.2 Stress Description
For a given layer s (= 1,2), made of an incompressible, isotropic and hyperelastic material, the
Cauchy stress tensor can be represented in polar coordinates as
θ eθ ⊗ eθ + T(s)

(2.22)
If the strain energy density is given as W(s) = W(s)(λr, λθ, λz), then the components of the

r er ⊗ er + T(s)

T(s) = T(s)

z ez ⊗ ez.

Cauchy stress tensor can be written in the following manner
∂W(s)
∂λq

T(s)
q = −π

(s) + λq

where π(s) are the Lagrangian multipliers, for both layers, introduced by the incompressibility
constraint, and q (= r, θ, z) signifies the three coordinate directions. From Eq. (2.21) we know that
the principal stretches are independent of θ and z, it follows that the Lagrangian multiplier can only
be function of r, namely π(s) = π(s)(r). Following a similar procedure to what was shown in Sec.
2.1, we have the strain energy function

For the case of plain strain, we can introduce a single variable strain energy function W(s) = ˆW(s)(λ)
such that λz = 1, λθ = λ and λr = 1/λ as follows

(2.23)

(2.24)

(2.25)

(2.26)

(2.27)

(2.28)

W(s) = W(s)(cid:0)λr, λθ, λz(cid:1)
(cid:19)
ˆW(s)(λ) = W(s)(cid:18) 1

, λ,1

λ

T(1)

r

r(1)

i

= T(2)

r

r(2)

e

= 0,

(cid:16)

(cid:17)

(cid:16)

(cid:17)

T(1)

r

r(1)

e

(cid:16)

(cid:17)

(cid:16)

(cid:17)

r(2)

i

.

= T(2)

r

26

In addition, for a cylindrical geometry, the equilibrium equations can be simplified as

dT(s)
r
dr

+

r − T(s)
T(s)

θ

r

= 0

which must be solved subject to the boundary conditions of zero traction at the inner and outer
faces,

as well as, the interface condition of perfect bonding between the two layers,

We employ the neo-Hookean strain energy function for modeling. This choice returns a form

for the Cauchy stress as

T(s) = −π

(s)
0 are the neo-Hookean coefficients for the sth layer. Therefore, we can write the components

here µ
of the stress tensor in the following manner

From the previous equation we get

T(s)

θ

= T(s)

(s)
r + µ
0

(cid:33)2(cid:35)

(2.30)

(2.31)

(2.29)

(2.32)

(2.33)

(2.34)

T(s)
r = −π

T(s)

θ

= −π

,

(s)
0 B(s)
(s)I + µ
(cid:33)2
(cid:32)
(cid:33)2
(cid:32)
(cid:33)2

R(s)
κr(s)
κr(s)
R(s)

(cid:34)(cid:32)

(cid:32)

,

(s)
(s) + µ
0

(s)
(s) + µ
0

κr(s)
R(s)

R(s)
κr(s)

−

(cid:32)

(cid:33)2

κr(s)
R(s)

(cid:34)(cid:32)
(cid:34)(cid:18) κr(s)

R(s)

−

R(s)
κr(s)

(cid:19)2

−

(cid:18) R(s)

κr(s)

(cid:33)2(cid:35)
(cid:19)2(cid:35)

dr(s)

∂T(s)
r
∂r(s) =

(s)
µ
0
r

eÞ

r

(s)

(s)
r
i

(s)
µ
0
r(s)

T(s)
r = κ

Like mentioned previously, since we are assuming a plain strain scenario, the stretch in the z
direction, λz = 1. Substituting these components of stress in the equilibrium equation, Eq. (2.26),
we obtain

which can be solved to give a value of radial Cauchy stress equal to

We can now rearrange Eq. (2.19) to substitute the value of R(s) into the previous equation. On
integrating (w.r.t boundary conditions in Eq. (2.27)), and simplifying, for our two layered problem
we get

(cid:34)
(cid:34)

(1)
T(1)
r = − µ
0
κ
(2)
T(2)
r = − µ
0
κ

2 log

κ

2 log

κ

(cid:33)
(cid:33)

(cid:32) R(1)
(cid:32)

i
R(1)
R(2)
e
R(2)

+ log

+ log

(cid:33)
(cid:33)

(cid:32)
(cid:32)

r(1)
r(1)
i
r(2)
r(2)
e

27

(cid:33)2(cid:41)(cid:35)
(cid:33)2(cid:41)(cid:35)

(cid:33)2
(cid:33)2

i

(cid:40)(cid:32) R(1)
(cid:40)(cid:32)

r(1)
R(2)
r(2)

e

i

e

(cid:32)
(cid:32)

−

−

R(1)
r(1)
R(2)
r(2)

+

+

1
2κ

1
2κ

θ

r

can be evaluated using Eq.

The corresponding values of T(s)
(2.31). Finally, we employed
these values of T(s)
in the interface condition, Eq. (2.28), to evaluate the radius in the deformed
configuration of the surface corresponding to the innermost surface in the reference configuration
for our system, r(1)
. The rest of the radii in the deformed configuration can subsequently be obtained
by using Eq. (2.20), and the fact that r(s+1)
2.3

Incremental Bifurcations

= r(s)
e .

i

i

Once the solution for the finite deformation is known, it can be used as the fundamental solution
for the analysis of incremental bifurcations. In the following we will describe the general incremen-
tal analysis applicable for any system in cylindrical coordinates, which is therefore applicable for
all cases of bending, unbending and eversion. We can define the incremental equilibrium equation

div(Σ) = 0,

where Σ is the incremental first Piola-Kirchhoff stress, defined as
, (cid:219)S = (cid:219)TF-T − TL-TF-T

Σ = (cid:219)SFT

.

(2.35)

(2.36)

Here, S = TF-T is the first Piola-Kirchhoff stress, and a superimposed dot indicates a first order
increment. In addition, the linearized constitutive equation is

Σ = CL − (cid:219)πI.

(2.37)

Here, C is the fourth-order tensor of instantaneous elastic moduli. In cylindrical coordinates, the
gradient of incremental displacements L can be written in the following manner

L = ur,r er ⊗ er +

eθ ⊗ eθ,
such that, when L is subject to the incompressibility constraint, i.e. tr L = 0, we get

er ⊗ eθ + uθ,r eθ ⊗ er +

uθ,θ + ur

r

r

ur,θ − uθ

rur,r + ur + uθ,θ = 0.

From Eq. (2.36) and the balance of rotational momentum we get

Σ12 − Σ21 = T2L12 − T1L21,

28

(2.38)

(2.39)

(2.40)

and on comparison with Eq. (2.37) we get (no summation over indices i and j)

Ci j ji + Ti = Cji ji,

(i (cid:44) j).

(2.41)

For a hyperelastic material, the components of C can be defined in terms of the strain energy
function, W. For our case of incompressible isotropic elastic materials, C can be written as a
function of two incremental moduli, denoted by µ and µ∗, that depends on the components of the
deformation gradient F. The non-zero components of C for the cylindrical coordinate system are

where

Crrrr = Cθθθθ = 2µ∗ + p,

Crθrθ = µ + Γ,

Cθrθr = µ − Γ,
Crθθr = Cθrrθ = µ + p,

Tθ − Tr

2

Γ =

and

p = −Tθ + Tr

2

,

(2.42)

(2.43)

describe the state of prestress. For a hyperelastic material, the two incremental moduli µ and µ∗
are given in terms of the strain energy function ˆW as

(cid:32)

(cid:33)

µ =

λ
2

4 + 1
4 − 1

λ

λ

d ˆW
dλ

,

µ∗ =

λ
4

(cid:32)d ˆW

dλ

(cid:33)

d2 ˆW
2
dλ

+ λ

(2.44)

(2.45)

For both scenarios of unbending and bending, we use the neo-Hookean material, for which the
incremental moduli µ and µ∗ are equal and simplify to

(cid:0)λ

4 + 1(cid:1)

2

λ

µ = µ∗ =

µ0
2

where the value of λ for bending is as given in Eq. (2.4), and for unbending and eversion it is given
as Eq. (2.21).

We seek bifurcations that can be represented by an incremental displacement field in the form

ur(r, θ) = f(r) cos nθ,
uθ(r, θ) = g(r) sin nθ,
(cid:219)π(r, θ) = k(r) cos nθ,

such that the incompressibility constraint given by Eq. (2.39) can be reformulated as

g = − f + r f (cid:48)

n

.

29

(2.46)

(2.47)

The substitution of the variable separated forms of the incremental displacement field, Eq.

(2.46), in the incremental equilibrium equation (div Σ = 0) gives the following results

(cid:18)
k(cid:48) = D f (cid:48)(cid:48) +
r2C
n2 f (cid:48)(cid:48)(cid:48) +

k =

(cid:19)
(cid:18) F
f (cid:48) +
n2 − D

(cid:19)

C + 2D

C,r + D,r +
F + 3C

r
r f (cid:48)(cid:48) +

n2

f ,

r2

E(1 − n2)
f (cid:48) − 1 − n2
n2

F
r

f ,

(2.48)

(2.49)

where, for a neo-Hookean material the coefficients C, D, E and F, are expressed as

C = µ − Γ =

µ0
2 ,
λ

D = 2µ∗ − µ =

E = µ + Γ = µ0λ

2

, F = rC,r + C

µ0
2

λ

4 + 1
λ

2

,

Note that for both bending and unbending, the coefficients C, D and E are equal, but the coefficient
F is different and given as follows
F = − µ0
2 ,
λ

(2.50)

(bending)
(unbending)

(cid:18)2

κ

(cid:19)

,

− 1
2

λ

F = µ0

No differentiation occurs while calculating C, D and E, which for us depend only on the type
of strain energy function employed, which is why they are equal. However, F is evaluated by
differentiating C w.r.t. r, which depends on the expressions for λ. Since λ is different for the finite
bending as compared to unbending, the expressions for F change accordingly.

By differentiating Eq. (2.48)2 with respect to r, and substituting the result in Eq. (2.48)1, we

get the following expression for f(r) within each layer

Cr4 f (cid:48)(cid:48)(cid:48)(cid:48)+2(F+2C)r3 f (cid:48)(cid:48)(cid:48)−(cid:0)(rF),r +4F−2n2D(cid:1)r2 f (cid:48)(cid:48)+(cid:0)(rF−2rn2D),r−2F(cid:1)r f (cid:48)+(1−n2)(F−rF,r−n2E) f = 0.

(2.51)
Once f(r) is known for each layer, the other functions g(r) and k(r) can be calculated by employing
Eqs. (2.47) and (2.48).

To find the set of all functions f (s)(r), with s = 1,2, we have to impose the following interface

and boundary conditions at the interfaces.

30

• Continuity of incremental tractions and displacements across interfaces

(cid:12)(cid:12)(cid:12)(cid:12)r=r
(cid:12)(cid:12)(cid:12)(cid:12)r=r

(1)
e

(1)
e

u(1)

r

(1)
rr

Σ

(cid:12)(cid:12)(cid:12)(cid:12)r=r
(cid:12)(cid:12)(cid:12)(cid:12)r=r

(2)
i

(2)
i

= u(2)

r

(2)
rr

= Σ

,

u(1)

θ

,

Σ

(1)
θr

(cid:12)(cid:12)(cid:12)(cid:12)r=r
(cid:12)(cid:12)(cid:12)(cid:12)r=r

= u(2)

θ

(2)
θr

= Σ

.

(2)
i

• Zero tractions at the boundaries defined by (r = r(1)

and r = r(2)
e )

i

(1)
e

(cid:12)(cid:12)(cid:12)(cid:12)r=r
(cid:12)(cid:12)(cid:12)(cid:12)r=r
(cid:12)(cid:12)(cid:12)(cid:12)r=r

(1)
e

(1),(2)
θr

(cid:12)(cid:12)(cid:12)(cid:12)r=r

(1),(2)
rr

Σ

= 0,

Σ

(2)
e

, r

(1)
i

= 0.

(2)
e

, r

(2.53)

(1)
i

(2)
i

,

(2.52)

• Shear stresses and incremental normal displacements vanish at the boundaries given by
θ = ± ¯θ, which gives

n =

2mπ
αl0

(m ∈ N).

(2.54)

Substituting the values from Eqs. (2.46), (2.47) and (2.48) for the variable separable equations,

in the equations for interface conditions, we get

(cid:110)

g

= g

(1)
e
(1)
e

(2)
i
(2)
i

f (1)(cid:12)(cid:12)(cid:12)r=r
= f (2)(cid:12)(cid:12)(cid:12)r=r
(1)(cid:12)(cid:12)(cid:12)r=r
(2)(cid:12)(cid:12)(cid:12)r=r
(cid:21)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)r=r
F − n2(cid:0)2D + C−Tr(cid:1)(cid:111)
F − n2(cid:0)2D + C − Tr(cid:1)(cid:111)
(cid:110)
(cid:21)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)r=r

r f (cid:48) + (n2 − 1)F f

=

(1)
e

Cr3 f (cid:48)(cid:48)(cid:48) + (F + 3C)r2 f (cid:48)(cid:48)+

Cr2 f (cid:48)(cid:48) + (C + Tr)r f (cid:48) + (n2 − 1)(C − Tr) f

(cid:20)

Cr3 f (cid:48)(cid:48)(cid:48) + (F + 3C)r2 f (cid:48)(cid:48) +

(cid:20)

(cid:20)

=

(1)
e

r f (cid:48) + (n2 − 1)F f

(2)
i

(cid:21)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)r=r
(cid:21)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)r=r

(2)
i(2.55)

Cr2 f (cid:48)(cid:48) + (C + Tr)r f (cid:48) + (n2 − 1)(C − Tr) f

(cid:20)

31

(cid:104)

Doing the same for zero tractions at external surfaces we get

Cr3 f (cid:48)(cid:48)(cid:48) + (F + 3C)r2 f (cid:48)(cid:48) +(cid:2)F−n2(2D + C − Tr)(cid:3)r f (cid:48) + (n2 − 1)F f
(cid:105)(cid:12)(cid:12)(cid:12)(cid:12)r=r

Cr2 f (cid:48)(cid:48) + (C + Tr)r f (cid:48) + (n2 − 1)(C − Tr) f

(cid:104)

(1)
i

(2)
e

, r

(cid:105)(cid:12)(cid:12)(cid:12)(cid:12)r=r

(1)
i
= 0

= 0

(2)
e

, r

(2.56)

In the next section, this above formulation will be solved numerically for a special case of a

bilayered system, undergoing bending, unbending and eversion deformation.

2.3.1 Compound matrix method for a bilayer

The compound matrix method was first proposed by Backus and Gilbert (1967). It was initially
employed to solve problems of fluid mechanics (Backus and Gilbert, 1967; Ng and Reid, 1979a,b,
1985), but was eventually used by Haughton and Orr (1995) for solving incremental elasticity in
the eversion of incompressible elastic cylinders. Later it was used to investigate the instability
in a homogeneous block (Haughton, 1999; Dryburgh and Ogden, 1999; Destrade et al., 2010)
undergoing flexure, as well as for an elastic multilayer subject to finite bending (Roccabianca et al.,
2011). Our approach to the compound matrix method is based heavily off the work of Roccabianca
et al. (2011) with only a few changes due to the notation and change in geometry of the problem.
The differential equation in Eq. (2.51) can be rewritten as a linear system of first order ODE’s.

The standard form for an elastic bilayer can be written in the following form

y(1)(cid:48)

= A(1)y(1)

,

y(2)(cid:48)

= A(2)y(2)

.

(2.57)

Here, the vectors y(s)(s = 1,2) , as well as the matrices A(s)(s = 1,2) are all functions of the
radial coordinate r, and the superscript (s) indicates the two layers in our system. The vectors
y(s)(s = 1,2) are defined as

f (s)(r)

f (s)(cid:48)(r)

f (s)(cid:48)(cid:48)(r)

32

(cid:104)

y(s)(r) =

f (s)(cid:48)(cid:48)(cid:48)(r)(cid:105)

(2.58)



The matrices A(s) are defined to be the following

A(s)(r) =

0
1
0
0
0
0
1
0
0
0
0
1
43 A(s)
42 A(s)
41 A(s)
A(s)
The components of A(s)(r) are given as follows
A(s)
41(r) =
A(s)
43(r) =
The boundary conditions on the external surface of the layer, with no loading (Eq. (2.56)) can

(n2 − 1)(F(s) − rF(s)
C(s)r4
2n2D(s) − (rF(s)),r − 4F(s)

2F(s) − (rF(s) − 2rn2D(s)),r

,r − n2E(s))

C(s)r2

(2.59)

(2.60)

44

,

,

A(s)
42(r) =
C(s)r3
44(r) = −2 F(s) + 2C(s)
A(s)
B(2)y(2)(cid:16)

C(s)r

r(2)

= 0,

e

be written in this form

B(1)y(1)(cid:16)

(cid:17)

r(1)

i

= 0,

B(s)(r) =

B(s)
(cid:1) = 0
(cid:1) + Hy(2)(cid:0)r(2)


,

where the matrices B(s)(s = 1,2) are written in the form

11 B(s)
B(s)
21 B(s)

12 B(s)
22 B(s)

23

13 B(s)
14
0
(cid:16)

Similarly, the interface equations can be written in a similar format, as follows
r(1)
e = r(2)

where

Gy(1)(cid:0)r(1)

e

,

i

(cid:17)

i
The matrices G and H are defined as follows

G =

G11 G12 G13 G14
0
0
0

G21 G22 G23
0
0

G41 G42

G31

0

H =

H11 H12 H13 H14
0
0
0

H21 H22 H23
0
0

H41 H42

H31

0

11 =(cid:0)n2 − 1(cid:1)F(s)
(cid:17)3
(cid:16)

The components of B(s)(s = 1,2) can be written as
B(s)
B(s)
14 =

B(s)
12 = r(s)
B(s)
21 = n2 − 1

F(s) − n2(cid:0)2D(s) + C(s)(cid:1)(cid:105)

B(s)
22 = r(s)

C(s)

r(s)

(cid:104)

o

o

o

(cid:16)
(cid:16)

B(s)
13 =
B(s)
23 =

o

r(s)
r(s)

o





33



(cid:17)


(2.61)

(2.62)

(2.63)

.

(2.64)

(cid:17)2(cid:0)F(s) + 3C(s)(cid:1)
(cid:17)2

(2.65)



,

1

r

r

i

,

,

(cid:16)
(cid:16)

(cid:17) are

(2.66)

G12 = rm

Here we define, r(1)

o = r(1)

, and r(2)

rm = r(i)

o = r(2)
e .
(cid:104)
(cid:104)

Similarly, the components of the matrices G and H such that(cid:16)
2 = r(e)
F(1) − n2(cid:16)2D(1) + C(1) − T(1)
(cid:17)(cid:105)
G11 =(cid:0)n2 − 1(cid:1)F(1)
n2(cid:16)2D(2) + C(2) − T(2)
(cid:17) − F(2)(cid:105)
H11 =(cid:0)1 − n2(cid:1)F(2)
G21 =(cid:0)n2 − 1(cid:1)(cid:16)
F(1) + 3C(1)(cid:17)
H21 =(cid:0)1 − n2(cid:1)(cid:16)
(cid:16)
F(2) + 3C(2)(cid:17)
(cid:17)
(cid:17)
(cid:16)

H12 = rm
G14 = r3
, H14 = −r3
G23 = r2
H23 = −r2

C(1) + T(1)
C(1) + T(1)

mC(1)
mC(2)
mC(1)
mC(1)

C(1) − T(1)
C(2) − T(2)

G13 = r2
m
H13 = −r2
m

,

,

,

r

r

G31 = 1,

G42 = rm
H42 = −rm

G41 = 1,
, H31 = −1, H41 = −1,

G22 = rm
H22 = −rm
Now, for each of the equations in Eq. (2.57) we have a family of solutions which satisfy each
equation. For each equation we can identify two independent solutions that satisfy each equations,
represented by y(s),I and y(s),II (s = 1,2). Due to how these equations were formulated, these
solutions satisfy the boundary conditions Eq. (2.61) but not the interface conditions Eq. (2.63).
For each layer, the two independent solutions of Eq. (2.57) can be arranged into two matrices, such
that their components sharing the following structure

(cid:17)
(cid:17)

r

r

,

,

which are defined as compound matrices for each layer s = 1,2. Now we introduce the vectors
(s)
k (s = 1,2, k = 1, . . . ,6), which are obtained by collecting the components of the minors of
φ
matrices in the following manner

(s),I
(s)
1 = y
1
(s)
(s),I
2 = y
1
(s),I
(s)
3 = y
1

φ

φ

φ

(s),II
y
2
(s),II
y
3
(s),II
y
4

(s),I
− y
2
(s),I
− y
3
(s),I
− y
4

(s),II
y
1
(s),II
y
1
(s),II
y
1

,

,

,

34

(s),I
(s)
(s),II
4 = y
y
2
3
(s)
(s),I
(s),II
5 = y
y
2
4
(s),II
(s),I
(s)
6 = y
y
4
3

,

(s),I
(s),II
− y
y
2
3
(s),I
(s),II
− y
,
y
4
2
(s),II
(s),I
− y
y
4
3

(2.68)

φ

φ

φ



(s),I
y
1
(s),I
y
2
(s),I
y
3
(s),I
y
4

(s),II
y
1
(s),II
y
2
(s),II
y
3
(s),II
y
4



(2.67)

From the work done in Ng and Reid (1979a) using Eqs.

(2.67) and (2.68), Eq.

(2.57) can be

rewritten in the following manner(cid:16)

φ(s)(cid:17)(cid:48)

= P(s)φ(s)

(s = 1,2)



0
0
0
0
1
P(s)
3

0
1
0
0
P(s)
2
−P(s)
1

0
0
1
1
P(s)
3
0



where P(s) is defined as the following matrix

0
1
P(s)
3
0
0
0
where the components of P(s) are given as follows

1
0
P(s)
2
0
0
P(s)
4

0
0
P(s)
1
0
P(s)
4
0

P(s) =

2F(s) −(cid:0)rF(s) + 2rn2D(s)(cid:1)
−2(cid:0)F(s) + 2C(s)(cid:1)

C(s)r3

C(s)r

,

P(s)
1 =
P(s)
3 =

2n2D(s) −(cid:0)rF(s)(cid:1)
(cid:0)1 − n2)(cid:1)(cid:0)F(s) − rF(s)

C(s)r2

,r − 4F(s)

,r − n2E(s)(cid:1)

,

C(s)r4

,r

,

P(s)
2 =
P(s)
4 =

(2.69)

(2.70)

(2.71)

(2.73)

The system of differential equations in Eq. (2.69) can be solved using a Runge-Kutta numerical
method to determine the compound matrices φ(s)(s = 1,2). The solution of this bifurcation problem
can be written as a linear combination of the solutions y(s),I and y(s),II,

y(1) = ξ1y(1),I + ξ2y(1),II,

y(2) = ξ3y(2),I + ξ4y(2),II

(2.72)

where the coefficients ξj (j = 1,2,3,4) are arbitrary, and represent the amplitude of the bifurcation
mode, and for this reason remain undefined in a linearized analysis. The conditions at the internal
interface Eq. (2.63) can be reformulated as

with ξ =(cid:2)ξ1 ξ2 ξ3 ξ4(cid:3)

Mξ = 0,

35

where

M =

(cid:16)Gy(1),I(cid:17)
(cid:16)Gy(1),I(cid:17)
(cid:16)Gy(1),I(cid:17)
(cid:16)Gy(1),I(cid:17)

1

2

3

4



(cid:16)Gy(1),II(cid:17)
(cid:16)Gy(1),II(cid:17)
(cid:16)Gy(1),II(cid:17)
(cid:16)Gy(1),II(cid:17)

(cid:16)Hy(2),I(cid:17)
(cid:16)Hy(2),I(cid:17)
(cid:16)Hy(2),I(cid:17)
(cid:16)Hy(2),I(cid:17)

1

2

3

4

1

2

3

4

(cid:16)Hy(2),II(cid:17)
(cid:16)Hy(2),II(cid:17)
(cid:16)Hy(2),II(cid:17)
(cid:16)Hy(2),II(cid:17)



1

2

3

4

such that, the bifurcation condition becomes

det(M) = 0

(2.74)

(2.75)

Since the values of the components of M depend on the various geometrical and material parameters
of the system there always exists a certain value of deformation for which the condition in Eq. (2.75)
is satisfied.

The previous equation can now be rewritten as the sum of 2 × 2-determinants as

1


i=0

(−1)i

(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)M1+i,1 M1+i,2

M42

M41

(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)

(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)M2−i,3 M2−i,4

M34

M33

−

M32

M12

M31

M14

M13

M43

M11

M44

(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)−

(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)M2+i,1 M2+i,2
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)M3−i,3 M3−i,4
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)M4−2i,3 M4−2i,4
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)M2+2i,1 M2+2i,2
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) = 0
2 + (Gk1Gl4 − Gk4Gl1(cid:1) φ
1 + (Gk1Gl3 − Gk3Gl1(cid:1) φ
5 + (Gk3Gl4 − Gk4Gl3(cid:1) φ
4 + (Gk2Gl4 − Gk2Gl4(cid:1) φ
2 + (Hk1Hl4 − Hk4Hl1(cid:1) φ
1 + (Hk1Hl3 − Hk3Hl1(cid:1) φ
5 + (Hk3Hl4 − Hk4Hl3(cid:1) φ
4 + (Hk2Hl4 − Hk2Hl4(cid:1) φ

(1)

(1)

(1)

(1)

(2.76)

(1)
3
(1)
6
(2.77)

These determinants can be expressed in terms of the compound matrices φ(s)(s = 1,2) as

Ml1 Ml2

(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Mk,1 Mk,2
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) =
and(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)Mk,3 Mk,4
(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12)(cid:12) =

(cid:0)Gk1Gl2 − Gk2Gl1(cid:1) φ
+(Gk2Gl3 − Gk3Gl2(cid:1) φ
(cid:0)Hk1Hl2 − Hk2Hl1(cid:1) φ
+(Hk2Hl3 − Hk3Hl2(cid:1) φ

(2)

Ml3 Ml4

(2)
3
(1)
6
(2.78)
where the indices k and l change values corresponding to their representation in Eq. (2.76). As an
example, we can see that setting k = 1 + i(i = 1,2) and l = 4 in Eq. (2.77) will give us the first half

(2)

(2)

(2)

36

of the first term in the summation of Eq. (2.76). Similarly, for the second half of the first term we
can set k = 2 − i(i = 1,2) and l = 3, and so forth

Once the geometrical parameters and material parameters of the system are fixed, the bifurcation
conditions (Eq. (2.75)), through the representation in Eq. (2.76), becomes a function of only the
deformation in the system, to be solved numerically. Due to how both the bending and unbending
deformations are described, for a given reference configuration, all deformations can be defined as
a function of the half angle of deformation. This is given by ¯θ for finite bending (see Fig. 2.1)
and αd for unbending and eversion (see Fig. 2.2). Therefore, for all cases this is the parameter for
which the instability problem is ultimately being solved for.

37

CHAPTER 3

USING FINITE BENDING TO EVALUATE PROPERTIES OF SINC

3.1 Methods

Here we take a novel approach to evaluate the mechanical properties of thin films of nanocrystals
on an elastomeric substrate. A PDMS-SiNC bilayer was used as a model to present the framework.
First, we inertially impacted a uniform film of SiNCs on PDMS. Then, we applied a finite bending
deformation to the system until instabilities formed on the surface coated with the SiNCs, noting
the critical angle at which the instabilities appeared. Finally, we applied a theoretical modeling
framework to describe the bilayer system. See Fig. 3.1 for an overview of our process.

3.1.1 PDMS Fabrication

To manufacture the PDMS, a Dow Corning heat-cured 184 Sylgard Elastomer PDMS kit was used
with a 10:1 weight ratio of base to curing agent. This mixture was stirred vigorously until air
bubbles formed, and placed under vacuum for 15 minutes to remove the trapped air. The mixture
was then weighed and poured into a petri dish of 85 mm diameter and was cured for 120 minutes at
60oC by placing on a hot plate to form PDMS of thickness 2.5 mm. After curing, the PDMS was
cut in rectangular sections according to our requirements and then attached to two metallic handles
with Gorilla® Super Glue as shown in Fig.3.1b and d.

3.1.2 Plasma synthesis of SiNCs

We produced the SiNCs using a nonthermal plasma reactor as has been reported previously (Mandal
and Anthony, 2016; Jariwala et al., 2011; Mangolini et al., 2005; Jurbergs et al., 2006). The reactor
consists of a pyrex tube through which we flow 94 sccm (standard cubic centimeters per minute) of
Argon (Ar) and 16 sccm of silane (SiH4, 5% in Ar), as well as 50 sccm of H2 through the sidearm for

38

Figure 3.1: Steps from plasma deposition to finite bending. SiNCs are deposited through a slit-
shaped orifice onto PDMS, which is then bent to observe formation of wrinkles on the SiNC
surface. (a) Photograph of plasma reactor. (b) Schematic diagram of reactor and SiNC synthesis.
(c) SiNC/PDMS sample illuminated with UV light, demonstrating photoluminescence. (d) Finite
bending applied with a lab-designed apparatus.

surface passivation. Radiofrequency (RF) power at 13.56 MHz was supplied to dual ring electrodes
encircling the reactor, and the pressure was kept constant at 3.85 Torr using a slit-shaped orifice.
The standoff distance between the orifice and the substrate was ∼4 mm. The downstream pressure
was ∼0.5 Torr. The nanocrystals were deposited in thin films by inertial impaction directly out
of the reactor onto the substrate(Anthony et al., 2012; Holman and Kortshagen, 2010; Nava et al.,
2018; Firth and Holman, 2018). One of the assets of this synthesis method is the room-temperature
process and ability to directly deposit SiNCs without post-processing steps or solvents, even on
temperature-sensitive substrates such as PDMS.

39

The PDMS samples were made (as described in Sec. 3.1.1), attached to a pushrod, and placed
under the orifice downstream of the reactor. The SiNCs were deposited by rastering the substrate
back and forth beneath the orifice to achieve a uniform film of SiNCs on the PDMS. Fig.3.1b shows
a work-flow schematic of the reactor and deposition of SiNCs on the PDMS. The rastering duration
was 20 minutes with each back and forth movement lasting 3 seconds, resulting in a SiNC layer
thickness of ∼4.5 µm.

Figure 3.2: SEM and FIB cross-sectioning of the SiNC layer. The thickness is estimated to be
(∼4.50 µm). (a) SEM image showing the surface of the SiNC layer on PDMS, including the FIB-
milled area (inset). (b) Higher magnification SEM image of the FIB-milled area demonstrating the
thickness of the SiNC layer.

40

Figure 3.3: (a) X-ray diffraction confirmed formation of SiNCs. XRD for SiNCs shows crystal
reflections at 2θ values 28.3, 47.3 and 56.1 degrees consistent with the (111), (220), and (311)
planes of silicon.
(b) TEM image confirming the formation of SiNCs: black circles highlight
individual SiNCs. The inset shows FFT of the correspondig SiNCs of the TEM image.

3.1.3 Silicon Nanocrystal Film Characterization

The thickness of the sample was measured by scanning electron microscope (SEM) imaging of
a cross-section of the sample generated using a focused ion beam (FIB). Fig.3.2a shows an SEM
image of the sample, and Fig.3.2b shows the FIB-generated cross section for thickness estimation.
To decrease charging effects and improve the image quality of the SiNC film, a thin layer of
Platinum (Pt) (∼6 nm) was deposited via Ar plasma sputtering prior to imaging. X-ray diffraction
(XRD), shown in Fig.3.3a, confirms the formation of crystalline nanoparticles. Transmission
electron microscopy (TEM) was performed on the SiNCs using a JEOL 2200 instrument, showing
individual NCs (Fig.3.3b). The TEM sample was prepared by tapping the TEM grid on the surface

41

of the sample to collect some of the SiNCs. The size of the SiNCs was measured using TEM image
analysis, giving an average SiNC diameter of ∼3.7 nm.

3.1.4 Uniaxial Tensile Test

We performed a standard uniaxial tensile test to evaluate the mechanical properties of the PDMS.
A rectangular sample of clear PDMS was mounted on a custom built uniaxial tensile test machine
and four black dots were applied on it as fiducial markers. The sample was then stretched cyclically
for 5 cycles from a configuration with a 10g pre-load (i.e., reference configuration) until a 20%
deformation was achieved. The sample was stretched between two clamps, one of which was
attached to a load cell to measure the force on line with the applied deformation. The strain rate
used during the tensile test was 0.1 mm/s. The force from the load cell was converted to engineering
stress by calculating the average cross-section area of the sample in the reference configuration. A
Hitachi KP-M2AN CCD camera was used to capture and feed the image of the fiducial markers
during uniaxial stretch into LabVIEW.

3.1.5 Finite Bending Experiment

We created a finite bending apparatus similar to the one used previously by Roccabianca et al.
(2010). Two important parameters are slenderness ratio (Λ = l0/h0) which is the ratio of width to
0 /hP0 ) which is the ratio of
thickness of the sample in cross section, and thickness ratio (H = hNC
thickness of the SiNC layer by that of the PDMS before finite bending (see Fig. 2.1a for reference).
The third dimension is taken to be 8-10 times larger than l0 to achieve a plane-strain condition on
the cross section. Here (and in all following equations) the superscript ‘P’ refers to the PDMS and
‘NC’ refers to the SiNCs.

The bilayer PDMS/SiNC system was put in the finite bending apparatus (see Fig.3.4) and
progressively bent until the SiNC film underwent bifurcation, at which point the critical angle was
recorded. Due to higher critical angles of bifurcation observed for samples with higher slenderness
ratios, deformation was applied by hand instead of using the apparatus. The onset of bifurcations

42

Figure 3.4: PDMS sample in the finite bending apparatus. (a) Side view, showing sample (with
handles) in the finite bending apparatus. (b) Top view showing formation of bifurcations.

was detected by direct visual inspection, and the critical angle was measured using analysis of
photographs of the side view of the sample. Each measurement was performed three times to
ensure accurate and reproducible readings.

3.1.6 Modeling bifurcations in elastic layered structures

Roccabianca et al. (2010), formulated a theoretical framework to study incremental bifrucations
due to finite bending of an elastic incompressible multilayered block. This theory was developed
for an N-layered system for Mooney-Rivlin type of material. The authors also performed some
experiments in bilayered structures that confirmed their theoretical predictions. Here, we employ a
similar theoretical modeling framework for the bilayered system to estimate the ratio of mechanical
properties of SiNC thin film to that of PDMS.

3.1.6.1 Finite deformation

The finite bending deformation is prescribed such that the reference configuration is rectangular
and described in a Cartesian coordinate system, and the deformed configuration is a sector of ring
described in a cylindrical coordinate system (see Fig. 2.1). Briefly, the finite bending deformation
can be described, within both the PDMS and the SiNC layers, by the following deformation gradient

43

F (Roccabianca et al., 2010)

F =

er ⊗ e

0
1 +

l0
2 ¯θr

2 ¯θr
l0

eθ ⊗ e

2 + ez ⊗ e
0

0
3,

(3.1)

where l0 is the overall width of the sample and 2 ¯θ is the angular deformation of the sample. We
consider the volume of the system to be conserved and the two layers to be perfectly bonded.
Therefore, we can apply the interface condition r NC
are the outer radius
of the SiNCs layer and the inner radius of the PDMS, respectively. This allows us to write,

, where r NC

o = r P
i

and r P
i

o

(cid:0)hNC

0

, hNC, α(cid:1) = r P

(cid:0)hP0 , hP, α(cid:1),

i

r NC
o

(3.2)

where hm
0 and hm represent the thickness of each layer before and after bending and the superscript
(m = NC, P), indicate the different layers in the system and α = 2 ¯θ/l0. Here, NC corresponds
to the nanocrystal layer, whereas P corresponds to the PDMS layer. To predict all geometrical
characteristics of the deformed state, we consider all the layers to be equilibrated. Briefly, we
define the Cauchy stress in each layer as,

T = −πI + F

T

,

F

dW
dC

(3.3)

T

F and I are the right Cauchy-Green deformation tensor and the identity tensor,
where C = F
respectively, and W and π are the strain energy function and the Lagrange multiplier for each layer,
respectively. We have employed the neo-Hookean strain energy function to describe both layers,
2 (trC − 3). Note that we are defining two constitutive equations, one for the
defined as W = µ0
PDMS layer as a function of the material parameter µP0 , and one for the SiNC layer as a function
of the material parameter µNC

.

0

For the layers to be equilibrated, the equilibrium equation in absence of body forces is required
to be satisfied within each layer. Furthermore, boundary conditions of zero traction are applied at
the innermost and outermost surfaces, as well as an interface condition of perfect bonding between
the two layers.

44

3.1.6.2 Bifurcation Analysis

The loss of uniqueness of the plane-strain incremental boundary value problem has been investigated
for multi-layered structures in previously published work (Roccabianca et al., 2010). Here we have
adapted the formulation to fit our application for a bilayered system of PDMS and SiNCs.

Briefly, for each layer, we seek to solve the incremental counterpart of the equilibrium, div (Σ) =
0, where Σ = (cid:219)SF
T is the incremental first Piola-Kirchhoff stress tensor, and the first Piola-Kirchhoff
-T for incompressible materials. Here,
stress and the Cauchy stress tensors are related by S = T F
and in the following, we use a superimposed dot to represent incremental quantities. The linearized
constitutive equation for each layer is given by

Σ = − (cid:219)πI + CL

(3.4)

where C is the fourth-order tensor of instantaneous elastic moduli and L = grad u is the gradient
of incremental displacement u(x). We consider bifurcations to be represented within each layer
by an incremental displacement field in a separable-variables form, as described by Roccabianca
et al. (2010). For incompressible isotropic elastic materials, one can write C as a function of two
incremental moduli, namely µ and µ∗. For neo-Hookean materials subject to finite bending in
plane strain condition (Bigoni, 2012), the elastic moduli can be written as

(cid:104)(αr)2 + (αr)−2(cid:105)

∗ =

µ = µ

µ0
2

.

(3.5)

Finally, the incremental boundary value problem is completed by a set of boundary and interface
conditions. Specifically, zero normal and shear incremental stresses at the innermost and outermost
layer, zero incremental shear stress and incremental radial displacement at the boundary θ = ± ¯θ,
and continuity of incremental stresses and displacements across the interface. To calculate the
angle that satisfies the bifurcation condition, θcr, we employed the compound matrix method, as
described in a previously published study (Roccabianca et al., 2011).

45

Figure 3.5: Critical angles of bifurcation evaluated as a function of ratio of material properties
for a PDMS-SiNC bilayer sample. The ratio of the neo-Hookean material parameters is defined
0 /µP0 . All curves have been calculated for a thickness ratio H included in the range
as ¯µ0 = µNC
0.0013–0.0018. Each line represent a different value of slenderness ratio, Λ = 2, black solid line;
Λ = 3, black dashed line; Λ = 4, dark gray solid line; Λ = 5, dark gray dashed line; Λ = 6, light gray
solid line. The insert represents a schematic of the bilayer and defines the slenderness ratio.

3.2 Results and Discussion

Experimentally, the critical angles of bifurcation were measured for nine samples by performing
the finite bending experiment described in Sec. 3.1.5. The values, reported in Table 3.1, have been
obtained for three experimental groups, characterized by a slenderness ratio Λ = 2, 3 and 5, each
group formed by a total of three samples (see Fig. 3.5a).

Employing the model described in Sec. 3.1.6, we evaluated the critical angle of bifurcation for
the bilayered structure, shown in Fig.3.5 and Fig.3.6. Specifically, after fixing the overall geometry
of the sample, namely the thickness ratio H and the slenderness ratio Λ, we calculated the value
of the angle representing the onset of bifurcation for each value of the neo-Hookean material

46

Table 3.1: Critical angles of bifurcation for different slenderness ratios, expressed in degrees.
Angles are reported for Λ = 2, 3 and 5.

Slenderness Ratio, Λ Critical angle, θcr Material Parameters ratio, ¯µ0

2
3
5

47.33 ± 3.21
71.33 ± 1.53
130.67 ± 3.21

1.90 ± 0.13
1.85 ± 0.05
1.68 ± 0.04

Figure 3.6: Plotting points from the experiments on graph obtained by numerical calculations. The
0 /µP0 . Shown here are points
ratio of the neo-Hookean material parameters is defined as ¯µ0 = µNC
obtained for Λ = 2, 3 and 5. Insert shows an instability on the surface of a sample with Λ = 2.

0 /µP0 .

In both Figures, we considered ¯µ0 to be included
parameters ratio, defined as ¯µ0 = µNC
between 0 and 100, however a pilot study showed that the value of the critical angles plateaus for
values of ¯µ0 included between 100 and 105 (data not shown). Furthermore, the same pilot study
showed that the critical angle of bifurcation is independent of the thickness ratio H, if H is included
in the range 0.0013 – 0.0018. The experimental value measured in this study is included within
that range, specifically H = 0.0017 on average for the nine samples.

47

Fig.3.5 shows the variation of the critical angle of bifurcation for varying values of the slen-
derness ratio, specifically for Λ included in the range 2 to 6. For a fixed value of ¯µ0, which is
represented by a vertical line in the graph, an increase in the slenderness ratio results in a later
appearance of the bifurcation, namely a higher value for the critical angle. This trend is consistent
with previously published results (Roccabianca et al., 2010, 2011).

Fig.3.6 shows how we employed the onset of bifurcation to estimate the mechanical properties
of the SiNC layer. Briefly, for each value of slenderness ratio considered in the experiments, we
calculated the correspondent theoretical sets of bifurcation angles. Then, the critical angle measured
experimentally was employed to uniquely identify the corresponding ¯µ0. The experimental value of
the critical angle fixed the abscissa of the point in the θcr – ¯µ0 plane, represented in Fig. 3.6. Then,
we imposed that the experimental point belonged on the bifurcation threshold curve calculated for
the corresponding slenderness ratio, which in turns fixed the ordinate of the point which is ¯µ0. The
value estimated for each experimental group, namely for each value of slenderness ratio considered,
is reported in Table 3.1. A vertical dashed line in Fig. 3.6 represents the average value of ¯µ0 for
all the experimental data, which is 1.81 ± 0.12. We then estimated the value of µP0 = 190 kPa by
fitting the data collected from the uniaxial tensile test, as described in Sec. 3.1.4. This gave us a
0 = 345 ± 23 kPa. The
value for the neo-Hookean material coefficient for the SiNC layer of µNC
neo-Hookean material constant is analogous to the shear modulus for low values of applied strain.
This finding is a clear indicator that the mechanical properties of this materials system are
determined by the interactions between SiNCs within the NC layer and between the NC layer and
the PDMS. Direct measurements of the mechanical properties of individual silicon nanocrystals
via microscope-guided nanoindentation have yielded elastic moduli around 170 GPa (Gerberich
et al., 2003), several orders of magnitude higher than our measurement of the shear modulus for
the SiNC layer on PDMS (assuming incompressible behavior, shear modulus is ∼1/3 of the elastic
modulus). Others have produced elastic/Young’s modulus estimates for single SiNCs between
70-170 GPa (Sadeghian et al., 2009; Mook et al., 2007; Kilymis et al., 2018). This range of
values (∼20-60 GPa for shear modulus) is also consistent with mechanical measurements of thin

48

ribbons of crystalline silicon on PDMS (Khang et al., 2006). By contrast, our value of 380 kPa
for shear modulus is far lower, and is consistent with the higher porosity of our nanopowder
layer as compared to bulk silicon. For example, Gross and Saber-Samandari (2009) showed that
increasing porosity of powder hydroxyapatite coatings leads to decreasing elastic modulus - further,
regardless of porosity, the coatings exhibited values of elastic modulus more than 35x smaller than
the modulus of single-crystal hydroxyapatite (Gross and Saber-Samandari, 2009; Zamiri and De,
2011). Similarly, the Young’s modulus of aluminum foams is between 10−1 and 10−3 the value
of Young’s modulus of solid aluminum, depending on porosity (Andrews et al., 1999). Finally,
previous results have shown that increased porosity significantly decreases the Young’s modulus
of polymer thin films deposited on PDMS, when measured by employing mechanical instabilities
(Stafford et al., 2004). The porosity of our layers was difficult to ascertain with high accuracy.
However, a preliminary density measurement provided us with a value near the limit for close-
packed randomly oriented rigid spheres (volume fraction 0.64). While we did not probe the
effects of porosity in this proof-of-concept investigation, the porosity of the SiNC layer is tunable
by adjusting reactor properties (Holman and Kortshagen, 2010; Firth and Holman, 2018), and
future work will be devoted to determining the relationship between layer porosity and mechanical
properties. This versatile method is broadly applicable to extract the mechanical properties of
nanoparticle layers on elastomer substrates regardless of material, layer thickness, and density. Our
results also offer a unique opportunity in the future to investigate how the NC surface functionality
(organic ligand coating, inorganic shell, etc.) influences the mechanical behavior of the NC films.

49

CHAPTER 4

UNBENDING AND EVERSION OF A BILAYER : RESULTS

The goal of this chapter is to discuss the results of unbending and eversion deformation discussed
in Sec. 2.2, as well as discussing the phenomenon of bifurcation associated with that deformation.
In the following we will focus on results collected for the case in which the reference configuration
is a full cylinder (as opposed to a sector of a cylinder), in other words the reference half-angle for all
cases considered here is π. However, the mathematical method presented in this thesis to evaluate
the onset of instability is a general formulation which can be applied to any referential angle.

In this chapter we first discuss the stress state for a few systems undergoing unbending and
eversion deformation, and discuss how changing the deformation in a system can influence the
complex stress state within it. Subsequently, we discuss the onset of instability due to deformation
in the system, and discuss how varying the various geometric and material parameters of the system
can influence the onset of bifurcation. Here we have specifically investigated how having a thin
stiffer layer on the inside, or outside of a ring influences the onset of bifurcation. We solve the
bifurcation problem for a variety of stiffnesses to evaluate the critical angle of deformation and the
connected stretch for the surface undergoing bifurcation.

4.1 Stress distribution in bilayer

In this section we will discuss the solutions obtained for the Cauchy stress distribution for a
bilayered system undergoing finite unbending / eversion deformation, as shown in Eqs. (2.31) and
(2.34). These analytical expressions, shown here for the first time, allow us to determine the complex
stress state within the cross section of a bilayered plate subject to this deformation with relative
ease. For example, it is intuitive to understand that, when the system undergoes unbending/eversion
deformation, the outermost radius in the reference state undergoes compression, while the innermost
radius in the reference state undergoes tension. However, this deformation can generate complex
stress distribution withing the thickness of the system, that could result in the presence of none,

50

one, or two neutral axis. To illustrate our point, in this section we will consider three different
systems, that under specific loading conditions show one, two, or zero neutral axis.

Figure 4.1: System undergoing unbending and eversion deformation. R(1)
e =
4.5, µr = 4. Deformed configurations are for αd = π/2,−π/2. The neutral axis for both deformed
configurations are as marked

= 2, R(2)

= 4, R(2)

i

i

i

= 4 and R(2)

In the first example, represented in Fig. 4.1, the system is characterized by the radii in the
e = R(2)
= 2, R(1)
reference configuration been equal to R(1)
e = 4.5, and by the ratio of
i
(1)
(2)
0 /µ
neo-Hookean coefficients being equal to µr = µ
0 = 4. For this system, we considered two
deformed states, one each for unbending and eversion, characterized by αd = π/2 (Fig. 4.1b) and
−π/2 (Fig. 4.1c). The radial stress, Tr, is compressive (negative) throughout the thickness, while
the circumferential stress, Tθ, can assume both tensile (positive) and compressive (negative) values
across the thickness. The location where the corumferential stress changes sign, in other words
assumes a value equal to zero, is named neutral axis. As expected, as the deformation increases,
moving from left to right in the Figure, the absolute peak values of both radial and circumferential
stresses become higher. Finally, the marked dashed line indicates the presence of a single neutral
axis for both loading conditions, as can be seen in Fig. 4.1 b and c.

In the second example, represented in Fig. 4.2, the system is characterized by the same

51

Figure 4.2: System undergoing unbending and eversion deformation. R(1)
e =
4.5, µr = 20. Deformed configurations are for αd = π/10,−π/2. The neutral axis for both deformed
configurations are as marked. Two neutral axes observed for deformation with αd = −π/2

= 2, R(2)

= 4, R(2)

i

i

i

i

= 2, R(1)

e = R(2)

= 4 and R(2)

geometrical parameters as in the previous case, namely the radii in the reference configuration are
taken to be equal to R(1)
e = 4.5, while the ratio of neo-Hookean
coefficients is taken to be equal to 20. For this system we have considered two loading conditions,
as described by values of αd = π/10 (bending), and αd = π/2 (eversion), as shown in Fig. 4.2
b and c. As before, in the Figure the marked dashed line indicates the presence of a neutral axis.
While in the case of unbending this system has one neutral axis, which happens for our case in the
inner layer, in the case of eversion we notice the presence of two neutral axes, one in each layer.
This is an interesting result that mimic what observed in the case of finite bending, as described
before in Roccabianca et al. (2010). As discussed in Sec. 1.2.1, although the presence of multiple

52

neutral axes had been reported before, (by Chu (1998) in thermal loading, and by Chuang and Lee
(2000) due to residual stresses) the result by Roccabianca et al. (2010) represented the first time
this had been reported in finite bending deformation. This here is the first report of multiple neutral
axes in the case of finite unbending / eversion deformation, to the best of our knowledge.

Figure 4.3: System undergoing unbending and eversion deformation. R(1)
e =
4.5, µr = 4. Deformed configurations are for αd = π/2,−π/2. The neutral axis for both deformed
configurations are as marked. Two neutral axes observed for deformation with αd = −π/2

= 3, R(2)

= 4, R(2)

i

i

In addition to the case in which two neutral axes were observed, we also found cases in which,
for a specific system in a specific loading condition, there were zero neutral axis. In other words,
in these conditions there were no locations across the thickness in which the stress along the
circumferential direction was equal to zero, intuitively this is only possible if one of the layers is
purely in compression whereas the other layer is purely in tension. One example of the case with
no neutral axis can be seen in Fig. 4.3. Specifically, we consider a system characterized by values
of radii in the reference configuration equal to R(1)
e = 4.5, and by the
value of the ratio of neo-Hookean coefficient taken equal to µr = 20. For this specific system, no
neutral axes are observed for an unbending deformation of αd = π/2, as shown in Fig. 4.3 b. For

= 4, and R(2)

e = R(2)

= 3, R(1)

i

i

53

the same system, a neutral axis can be observed for the deformation characterized by αd = −π/2,
as shown in Fig. 4.3 c.

The only results that we could find in literature which mentioned the absence of neutral axes
during bending was in the work by Hsueh (2002) and Hsueh et al. (2003). They discussed that
while solving for bending deformations in bilayers with residual stresses they were able to find
certain scenarios in which no neutral axes were observed. Therefore, the result that we obtain here
is a novel result not only because this is the first time that systems have been observed undergoing
unbending and eversion deformation with no neutral axes, but also this is the first time that we
observe no neutral axis for a system with no residual stresses.

4.2 Bifurcation behavior during unbending and eversion

4.2.1 Single layer

Figure 4.4: Critical deformation angles αd in unbending and eversion, plotted versus ratio of
innermost to outermost radius ratio ρ. Note that, this study was done for a single-layered system,
and that the reference configuration had varying angles.‘m’ refers to the mode number of bifurcation
(Sigaeva et al., 2018)

As discussed in Sec. 1.3 one previous study focused on the bifurcation analysis of the problem
of a single-layered sector of hollow circular cylinder undergoing bending, unbending and ever-

54

sion (Sigaeva et al., 2018). The authors studied the onset of bifurcation for varying reference
configurations, namely for varying innermost to outermost radius ratio (ρ), using the impedance
matrix method. Furthermore, they considered several reference angles (αr), to analyze bifurcation
behavior for a deformation angle (αd) from αr to −π. Since in this section we focus on the defor-
mations of unbending and eversion, we report a compilation of the results they obtained for these
deformations in Fig. 4.4 (results for bending not shown). The figure shows the envelope of the
bifurcation angles for seven different reference angles, varying from αr = π/12 to αr = π. For each
geometry, and for each value of ρ only the mode number (m) corresponding to the highest critical
angle of bifurcation is reported. Henceforth, this mode number is referred to as the primary mode
number. The reason we only note the primary mode number is because, for a system with given
parameters, theoretically, we can calculate bifurcations corresponding to multiple different mode
numbers. But, the only mode number which can be observed physically is the first mode number
encountered in deformation. Since, for our geometry the system deforms from an angle of π to −π
the mode number with the maximum value of critical angle of bifurcation is reached first, and is
the mode number of significance for that particular system.

Note that, the minimum value that αd can assume is −π. For this value of deformation the
system is completely everted, which means that it forms a tube in which the inner and outer faces
invert roles when compared to the reference configuration. For this reason, the lower limit of the
vertical axis in Fig. 4.4 is equal to −π. Finally, from the figure we can extrapolate that, for each
value of reference angle, there exist a limiting value of ρ above which the system will be able to
undergo complete eversion without bifurcations.

In this study, we have considered bilayered systems with an initial reference angle of αr = π,
namely the refernce configuration is always a complete hollow cylinder. First, we evaluated the
bifurcation behavior for a system for which both layers are made of the same material. That allows
us to replicate the condition of a single-layered system so that it can be compared to what provided
by Sigaeva et al. (2018), to verify the precision of our analysis. Fig. 4.5 represents a comparison
between the solution presented in this study, in the conditions just described, and the solution

55

Figure 4.5: Studying the variation of critical angle of bifurcation θcr with the ratio of innermost to
outermost ratio of the system ρ. Comparing results published by Sigaeva et al. (2018) for unbending
in a single layer, to our own results obtained by employing compound matrix method on a bilayer.
Results for the single layer were replicated by applying the same material properties to both layers

presented in Sigaeva et al. (2018) for a reference angle αr = π.

At low values of ρ, the solutions obtained by both methods is remarkably similar, as expected.
However, a difference between the two solutions is observed, at high values of ρ. This seems to be
associated with the re-surfacing of the first mode number at high values of ρ (slender structures),
as a primary mode of bifurcation (as shown in Fig. 4.6). This was not reported by Sigaeva et al.
(2018). A similar behavior, however, was discussed by Roccabianca et al. (2011) for the case
of bifurcation in bi-layered structures subject to bending. The authors reported that, for certain
combinations of initial geometries and material properties, the first mode number while not being
the primary mode at low slenderness ratio, could become the primary mode at higher slenderness
ratios. One possible explanation of the difference between the results presented here and the results
reported in Sigaeva et al. (2018) is that the two studies employ different numerical methods to solve

56

Figure 4.6: Solving the bifurcation problem to study how instability in different mode numbers
behaves as a function of innnermost to outermost ratio of the system ρ for a single-layered system
undergoing unbending. Solution for mode numbers from 1 to 12 is presented and the mode number
corresponding to highest critical angle of bifurcation is marked for all ρ to indicate the first mode
which will be observed in a real world deformation. Note the fact that mode number 1 is the
significant mode both at very small, and for higher values of ρ

the incremental equilibrium problem. Here, we employed the compound matrix method while the
authors of Sigaeva et al. (2018) used the impedance matrix method.

4.2.2 Bilayered systems

A non-exhaustive list of the various different geometrical parameters that can be associated with
our bilayered system in the reference configuration are as follows

• the inner and outer radii of the layers (R(1)

i

, R(2)

i

, R(1)
e ),

• the thickness of the individual layers,

• the ratio of radii and thicknesses of the layers.

57

However to uniquely identify the geometry in the reference configuration we need to define only
three appropriate parameters. For example, we can uniquely identify a geometry if we select the
three radii of the system, or if we choose two thicknesses corresponding to both layers and any one
radius of the system, or any different variety of combinations of the above. To study the variation
of angle of bifurcation with varying reference configurations, we studied a family of bilayers with
one layer of unit thickness and the other with a thickness of 0.2. The only other parameters we
need to uniquely define a system is the ratio of neo-Hookean coefficients, µr of the system, here
defined as the ratio of the coefficients for the outer layer to the inner layer, as seen in the reference
(1)
(2)
0 /µ
0 according to the terms defined in Sec. 2.2). We then varied the value
configuration (µr = µ
of the inner radius, which resulted in a variation of ρ.

The stiffness of the thicker layer is kept constant and the thinner layer is varied to observe
the difference in bifurcation behavior due to changing µr. Following that, the critical angle of
bifurcation is numerically evaluated as a function of the ratio of the innermost to outermost radius,
ρ. Finally, the critical stretch (for the layer corresponding to the outermost layer in reference
configuration) is also evaluated and reported.

First, we will discuss results calculated for a system characterized by the thinner layer placed
at the outermost radius and a ratio of material parameters µr = 10, as presented in Fig. 4.7. In
other words, we are considering a thick hollow tube made of a soft material, that is coated on the
outer surface by a thinner layer (20% of the original thickness) made of a stiffer material (10 times
the original stiffness). In the Figure, we represent the values of the critical angle of bifurcation and
of the circumferential stretch on the bifurcated surface (surface in compression) as a function of
the ratio of innermost to outermost radius, ρ. As observed in the case of µr = 1, the first mode
(m = 1) becomes critical at higher values of ρ, however, at lower values of ρ, the primary mode
corresponds to m = 3, as opposed to m = 1. As far as we know, no one has reported such a scenario
in which a higher mode number than 1 is the primary mode number at higher slenderness. Another
very interesting trend that can be observed is the fact that for all values of ρ the bilayer becomes
unstable earlier (i.e., at a higher critical angle and stretch of bifurcation) than the single-layered

58

Figure 4.7: Behavior of stretch and angle of deformation for a bilayered system as a function of
innnermost to outermost ratio of the system ρ. Inner layer of the system has unit thickness, outer
layer a thickness of 0.2 and ratio of inner to outer neo-Hookean coefficients is 10. Similar to
Fig.4.6, the mode number corresponding to highest critical angle of bifurcation is marked for all ρ
to indicate the first mode which will be observed in a real world deformation

system. Therefore, the addition of a stiff layer outside can drastically affect the bifurcation behavior
of the system. For all values of ρ the bifurcation happens at a value significantly higher than the
compressive stretch for surface instability in an elastic half space (λcr = 0.543). This is quite

59

similar to what was observed in the analysis of instability in bilayers, with a large difference in
neo-Hookean coefficients, undergoing finite bending in the case where the stiffer surface was the
one experiencing bifurcation.

In the next sections we will discuss the effect of a thin, stiff layer added at the compressive
(outermost) or tensile (innermost) sides of the thick hollow cylinder considered here. We will also
consider the effect of different stiffness ratios between the two layers. Numerical solutions are
provided for values of cases in which the stiff layer is 5, 10, 30, 50 and 80 times stiffer both when
it’s on the compressive and tensile side. In both cases considered, we included the case of µr = 1 to
replicate the results for the onset of bifurcation in a single layered system, for ease of comparison.

4.2.2.1 Thinner, stiffer layer at the compressed side (outermost radius)

In Fig.4.8, we show the results, in terms of critical circumferential stretch (top) and critical angle of
bifurcation (bottom), calculated for the case of a cylinder made of a layer of unit thickness, which has
been coated at the outer side with a thinner layer, with thickness 0.2. Like mentioned previously,
the thinner layer is also taken to be stiffer when compared to the thicker layer; specifically, we
considered different values of µr (= 1,5,10,30,50,80), to analyze the effect of a stiffening outer
thin layer on the bifurcation behavior of this system undergoing unbending and eversion.

Additionally, Fig. 4.8 shows that, with an increase in the stiffness of the thinner outer layer,
there is a general shift upwards of the value of critical circumferential stretch. This means that a
bilayered system becomes unstable at progressively lower values of deformation applied, resulting
in higher critical angle θcr as well as critical circumferential stretch λcr, as the outer layer becomes
progressively stiffer. This is especially true for low values of ρ where there is a direct upward trend
for θcr and λcr for bifurcation.

However, not all aspects of the results show a monotonic behavior, with increasing values of
µr. This can be seen in the behavior of the critical angle of bifurcations for certain values of
ρ. To show what we mean, let’s take an example. For ρ = 0.3, the critical angle of bifurcation
increases when µr goes from 1 to 10, but on further increasing µr to 30 and onwards the critical

60

(a)

(b)

Figure 4.8: Variation in bifurcation behavior for bilayered system with outer layer thickness = 0.2,
unit thickness for inner layer, and ratio of inner to outer layer material properties as indicated in
the legend. The variation was observed as a function of the innermost to outermost radius ratio (ρ)
for two parameters (a) the stretch, and (b) the angle of deformation, at which the system becomes
unstable. The dots indicate the point at which a different mode number becomes the primary mode
of bifurcation

61

angle of bifurcation seems to drop sharply. At even higher values of ρ (cid:39) 0.45 the critical angle of
bifurcation for a system with µr = 5 is higher than a system which has µr = 10. This seems to
suggest that, contrary to what we found for low values of ρ, where the critical angle of bifurcation
monotonically increased with increasing µr, for higher values of ρ we instead find that there exists
a particular value of µr for which we attain the highest critical angle of bifurcation. At any value
of µr higher than that particular value, the system again becomes unstable at a smaller angle.

On one hand, this analysis shows that the bifurcation behavior of hollow cylinder undergoing
unbending and eversion is significantly impacted by the application of a thin external coating made
of a stiffer material. On the other hand, there seems to be a upper limit on the value of ρ, above
which the system can be everted completed in a stable manner, and this limit seems to be somewhat
independent of the stiffness of the outer coating. For the geometry considered here, this limit is
ρ ≈ 0.6. This limit could exist because we chose a thickness ratio of outer to inner layer to be
0.2 for our system. It is possible that changing the thickness ratio could bring about a change in
this limit, however more data needs to be calculated to make any solid conclusions on this matter.
Additionally, based on our discussion of Sigaeva et al. (2018) it is reasonable to assume that varying
the reference angle of the system (αr) would bring about a change in this limit.

4.2.2.2 Thinner, stiffer layer at the tensile side (innermost radius)

In Fig. 4.9 we show the critical circumferential stretch (top) and angle of bifurcations (bottom)
evaluated for a hollow cylinder of unit thickness undergoing unbending/eversion, coated at the inner
radius with a thinner layer (thickness of 0.2) made of stiffer material. The material parameters
ratio between the two layers was varied, specifically µr (= 1,1/5,1/10,1/30,1/50,1/80). This
problem somewhat mimics what was considered in the previous example (Sec. 4.2.2.1), the main
difference being that here the stiff coating is considered to be at the tensile side (innermost radius
in the reference configuration) while in the previous section we considered the coating to be at the
compressed side (outermost radius in the reference configuration).

Fig. 4.9, also indicates that the addition of a progressively stiffer layer at the innermost radius

62

(a)

(b)

Figure 4.9: Variation in bifurcation behavior for bilayered system with inner layer thickness = 0.2,
unit thickness for outer layer, and ratio of inner to outer layer material properties as indicated in
the legend. The variation was observed as a function of the innermost to outermost radius ratio (ρ)
for two parameters (a) the stretch, and (b) the angle of deformation, at which the system becomes
unstable. The dots indicate the point at which a different mode number becomes the primary mode
of bifurcation

63

of a hollow cylinder results in a progressive increase of both λcr and θcr. In other words, for a
given value of innermost to outermost radius ratio ρ, the sample with stiffer inner layer becomes
unstable always for smaller value of applied deformation, resulting in a larger critical angle as
well as compressive circumferential stretch, when compared to a samples with softer inner layer
and to the sample made of one material. An interesting thing to note is that for values of ρ
approaching 0, all samples become unstable for approximately the same value of unbending angle
(θcr ≈ 0.9). Furthermore, for low values of µr (i.e., higher stiffness of thin, inner layer), the critical
circumferential stretches also seem to be converging to a constant when ρ approaches 0. They
actually seem to tend very close to the value associated with surface bifurcation λcr = 0.543. This
is actually consistent with similar results for bifurcation in a bilayer undergoing finite bending, as
reported by Roccabianca et al. (2010). For a case in which a stiff layer was on the tensile side,
the authors reported bifurcation behavior tending towards surface instability for low slenderness
geometries.

As the values of ρ approach 1, for high values of µr, λcr seems to keep increasing due to the
behavior of mode m = 1, until we reach a deformation of αd = −π. Interestingly, for higher values
of ρ the lines indicating bifurcation behavior seem to be converging to a limit. This trend can be
noticed in both the graphs for λcr and θcr vs ρ where at lower, the lines for different µr start getting
very close to one another. In other words, there is a smaller change in bifurcation behavior when
changing the stiffness of the inner layer from 1 to 10, than it is when changing it from 50 to 80.
This would suggest that there are diminishing returns, as far as change in bifurcation behavior that
can be observed by increasing the µr is concerned, for very high values of µr.

Finally, in the previous Section, we have learned that adding a coating on the outer side of the
cylinder did not affect the lower limit of ρ for which the systems can be completely everted in a
stable manner, independently of the value of the stiffness of the coating. For the specific thickness
ratio considered here, that limit value of ρ was approximately equal to 0.6. However, Fig. 4.9
shows that, in the case of adding a coating to the cylinder at the inner side, the limiting value of ρ
above which it is possible to evert the system in a stable manner depends on the value of stiffness

64

of the coating. Specifically, we notice that an increase in the stiffness of the coating added at the
inner side leads to an increase in the limit value of ρ. In other words, the stiffer the inner layer, the
more slender a structure is needed to be able to evert without loss of stability.

65

CHAPTER 5

FUTURE WORK

The method to evaluate the mechanical properties of thin films in SiNC as outlined in Chap. 3 is
a versatile technique and can be employed to evaluate the mechanical properties of any thin film
deposited or laid on a soft substrate. However, due to the novelty of this method it becomes essential
to try and develop another technique that can be used to evaluate mechanical properties for such a
system. The work presented in Sec. 2.2 along with the results presented in Chap. 4 provide a few
essential first steps towards formulating a method with which we can validate the results obtained
in Chap. 3.

5.1 Evaluating the mechanical properties of thin films

deposited on a pre-stretched surface

A discussion regarding the kinematics behind unbending a bilayered cylindrical system has
been dealt with in Sec. 2.2. The system discussed there is based on the assumption that the bilayer
just exists as a cylinder and unbending deformation is carried out to observe the onset on instability.
However the system we have in real life includes an additional step beforehand which involves
wrapping a rectangular PDMS sample around a cylinder. Therefore, the problem we need to work
on to evaluate mechanical properties in unbending, are slightly more complicated, and are detailed
in Fig. 5.1.

The thickness and width of the PDMS sample, in the reference configuration, is taken in such
a way that there are no tractions on the surface when it is rolled on a cylindrical rod of a given
outer radius. On this stress-free surface the SiNC are deposited on the PDMS substrate. F1 and
F2 are the deformation gradients obtained on applying the bending and unbending deformations,
respectively (see Fig. 5.1). F = F1F2 is the deformation gradient for the total deformation in the
system, from the reference configuration. In the special case in which there is no deposition on the

66

Figure 5.1: Schematic showing the process by which samples are made for evaluating material
properties of a SiNC deposited on a pre-stretched substrate. The substrate, which is initially
rectangular, is rolled into a cylinder. The SiNC is deposited on it using the plasma deposition
process outlined in Sec3.1. The sample is then opened up and bifurcations are observed on the
compressed surface of the deposited layer

surface, and the system unbends to F2 = F−1
deformation of the bending. In that case F = F1F−1
has taken place in the system.

1 , i.e., the unbending deformation F2 is just the inverse
1 = I, which implies that overall no deformation

When the system undergoes the unbending deformation the thin SiNC layer deposited on the
outside surface of the PDMS experiences a compressive force. This is similar to the compressive
force faced by the top surface of the SiNC in the method involving finite bending. Therefore,
like in the finite bending experiments, we expect the formation of instabilities on the surface of
thin film above a critical compressive stress. This value of compressive stress depends upon the
geometric configuration and mechanical properties of the system. We aim to use the formation
of these instabilities to evaluate and validate the mechanical properties of thin film SiNC when
deposited on a soft PDMS substrate.

67

BIBLIOGRAPHY

68

BIBLIOGRAPHY

Andrews, E., Sanders, W., and Gibson, L. J. (1999). Compressive and tensile behaviour of aluminum

foams. Materials Science and Engineering A, 270(2):113–124.

Anthony, R. J., Cheng, K.-Y., Holman, Z. C., Holmes, R. J., and Kortshagen, U. R. (2012). An
all-gas-phase approach for the fabrication of silicon nanocrystal light-emitting devices. Nano
letters, 12(6):2822–2825.

Anthony, R. J., Rowe, D. J., Stein, M., Yang, J., and Kortshagen, U. (2011). Routes to achieving
high quantum yield luminescence from gas-phase-produced silicon nanocrystals. Advanced
functional materials, 21(21):4042–4046.

Aron, M. (2000). Some remarks concerning a boundary-value problem in nonlinear elastostatics.

Journal of elasticity and the physical science of solids, 60(2):165–172.

Aron, M. (2005). Combined axial shearing, extension, and straightening of elastic annular cylin-

drical sectors. IMA journal of applied mathematics, 70(1):53–63.

Aron, M., Christopher, C., and Wang, Y. (1998). On the straightening of compressible, nonlinearly

elastic, annular cylindrical sectors. Mathematics and Mechanics of Solids, 3(2):131–145.

Arruda, E. M. and Boyce, M. C. (1993). A three-dimensional constitutive model for the large
stretch behavior of rubber elastic materials. Journal of the Mechanics and Physics of Solids,
41(2):389–412.

Backus, G. E. and Gilbert, J. (1967). Numerical applications of a formalism for geophysical inverse

problems. Geophysical Journal International, 13(1-3):247–276.

Bao, Z. and Chen, X. (2016). Flexible and Stretchable Devices. Advanced Materials, 28(22):4177–

4179.

Bigoni, D. (2012). Nonlinear solid mechanics: bifurcation theory and material instability. Cam-

bridge University Press.

Biot, M. A. (1963). Surface instability of rubber in compression. Applied Scientific Research,

Section A, 12(2):168–182.

Biot, M. A. (1964). Mechanics of incremental deformations.
Bowden, N., Brittain, S., Evans, A. G., Hutchinson, J. W., and Whitesides, G. M. (1998). Spon-
taneous formation of ordered structures in thin films of metals supported on an elastomeric
polymer. Nature, 393(6681):146.

Chu, S. (1998). Elastic bending of semiconductor wafer revisited and comments on stoney’s

equation. Journal of the Electrochemical Society, 145(10):3621–3627.

69

Chuang, T. J. and Lee, S. (2000). Elastic flexure of bilayered beams subject to strain differentials.

Journal of Materials Research, 15(12):2780–2788.

Chung, J. Y., Nolte, A. J., and Stafford, C. M. (2011). Surface wrinkling: a versatile platform for

measuring thin-film properties. Advanced Materials, 23(3):349–368.

Coman, C. D. and Destrade, M. (2008). Asymptotic results for bifurcations in pure bending of

rubber blocks. Quarterly Journal of Mechanics and Applied Mathematics, 61(3):395–414.

Destrade, M., Gilchrist, M. D., Motherway, J. A., and Murphy, J. G. (2010). Bimodular rubber

buckles early in bending. Mechanics of Materials, 42(4):469–476.

Destrade, M., Ogden, R., Sgura, I., and Vergori, L. (2014a). Straightening wrinkles. Journal of the

Mechanics and Physics of Solids, 65:1–11.

Destrade, M., Ogden, R. W., Sgura, I., and Vergori, L. (2014b). Straightening: existence, uniqueness

and stability. Proceedings of the Royal Society A, 470:20130709.

Dryburgh, G. and Ogden, R. W. (1999). Bifurcation of an elastic surface-coated incompressible
isotropic elastic block subject to bending. Zeitschrift fur Angewandte Mathematik und Physik,
50(5):822–838.

Ericksen, J. L. (1954). Deformations possible in every isotropic, incompressible, perfectly elastic

body. Zeitschrift für angewandte Mathematik und Physik (ZAMP), 5(6):466–489.

Fedorchenko, A. I., Wang, A.-B., and Cheng, H. H. (2009). Thickness dependence of nanofilm

elastic modulus. Applied Physics Letters, 94(15):152111.

Firth, P. and Holman, Z. C. (2018). Aerosol Impaction-Driven Assembly System for Production
of Uniform Nanoparticle Thin Films with Independently Tunable Thickness and Porosity. ACS
Applied Nano Materials, 1:acsanm.8b01334.

Gaspar, J., Paul, O., Chu, V., and Conde, J. (2010). Mechanical properties of thin silicon films
deposited at low temperatures by pecvd. Journal of Micromechanics and Microengineering,
20(3):035022.

Gei, M., Bigoni, D., and Guicciardi, S. (2004). Failure of silicon nitride under uniaxial compression

at high temperature. Mechanics of materials, 36(4):335–345.

Gent, A. (1996). A new constitutive relation for rubber. Rubber chemistry and technology,

69(1):59–61.

Gent, A. N. and Cho, I. S. (1999). Surface Instabilities in Compressed or Bent Rubber Blocks.

Rubber Chemistry and Technology, 72(2):253–262.

Gerberich, W. W., Mook, W. M., Perrey, C. R., Carter, C. B., Baskes, M. I., Mukherjee, R., Gidwani,
A., Heberlein, J., McMurry, P. H., and Girshick, S. L. (2003). Superhard silicon nanospheres.
Journal of the Mechanics and Physics of Solids, 51(6):979–992.

70

Greer, J. R., Oliver, W. C., and Nix, W. D. (2005). Size dependence of mechanical properties of
gold at the micron scale in the absence of strain gradients. Acta Materialia, 53(6):1821–1830.
Gross, K. A. and Saber-Samandari, S. (2009). Revealing mechanical properties of a suspen-
sion plasma sprayed coating with nanoindentation. Surface and Coatings Technology, 203(20-
21):2995–2999.

Han, X., Zheng, K., Zhang, Y., Zhang, X., Zhang, Z., and Wang, Z. L. (2007). Low-temperature in

situ large-strain plasticity of silicon nanowires. Advanced Materials, 19(16):2112–2118.

Haughton, D. (1999). Flexure and compression of incompressible elastic plates. International

journal of engineering science, 37(13):1693–1708.

Haughton, D. and Orr, A. (1995). On the eversion of incompressible elastic cylinders. International

journal of non-linear mechanics, 30(2):81–95.

Haughton, D. and Orr, A. (1997). On the eversion of compressible elastic cylinders. International

journal of solids and structures, 34(15):1893–1914.

Holman, Z. C. and Kortshagen, U. R. (2010). A flexible method for depositing dense nanocrystal

thin films: impaction of germanium nanocrystals. Nanotechnology, 21(33):335302.

Hsueh, C., Lee, S., and Chuang, T. J. (2003). An alternative method of solving multilayer bending
problems. TRANSACTIONS-AMERICAN SOCIETY OF MECHANICAL ENGINEERS JOUR-
NAL OF APPLIED MECHANICS, 70(1):151–153.

Hsueh, C.-H. (2002). Modeling of elastic deformation of multilayers due to residual stresses and

external bending. Journal of Applied physics, 91(12):9652–9656.

Jariwala, B. N., Dewey, O. S., Stradins, P., Ciobanu, C. V., and Agarwal, S. (2011). In situ gas-phase
hydrosilylation of plasma-synthesized silicon nanocrystals. ACS applied materials & interfaces,
3(8):3033–3041.

Jiang, C., Singamaneni, S., Merrick, E., and Tsukruk, V. V. (2006). Complex buckling insta-
bility patterns of nanomembranes with encapsulated gold nanoparticle arrays. Nano Letters,
6(10):2254–2259. PMID: 17034093.

Jurbergs, D., Rogojina, E., Mangolini, L., and Kortshagen, U. (2006). Silicon nanocrystals with

ensemble quantum yields exceeding 60%. Applied Physics Letters, 88(23):233116.

Khang, D. Y., Jiang, H., Huang, Y., and Rogers, J. A. (2006). A stretchable form of single-crystal

silicon for high-performance electronics on rubber subtrates. Science, 311(5758):208–212.

Kilymis, D., Gérard, C., Amodeo, J., Waghmare, U., and Pizzagalli, L. (2018). Uniaxial compres-
sion of silicon nanoparticles: An atomistic study on the shape and size effects. Acta Materialia,
158:155–166.

Lockwood, A. and Inkson, B. (2008). In situ tem nanoindentation and deformation of si-nanoparticle

clusters. Journal of Physics D: Applied Physics, 42(3):035410.

71

Lu, C., Dönch, I., Nolte, M., and Fery, A. (2006). Au nanoparticle-based multilayer ultrathin films
with covalently linked nanostructures- spraying layer-by-layer assembly and mechanical property
characterization. Chemistry of materials, 18(26):6204–6210.

Mandal, R. and Anthony, R. J. (2016). Aging of silicon nanocrystals on elastomer substrates:

Photoluminescence effects. ACS applied materials & interfaces, 8(51):35479–35484.

Mangolini, L., Thimsen, E., and Kortshagen, U. (2005). High-yield plasma synthesis of luminescent

silicon nanocrystals. Nano letters, 5(4):655–659.

Mook, W. M., Nowak, J. D., Perrey, C. R., Carter, C. B., Mukherjee, R., Girshick, S. L., McMurry,
P. H., and Gerberich, W. W. (2007). Compressive stress effects on nanoparticle modulus and
fracture. Physical Review B - Condensed Matter and Materials Physics, 75(21):1–10.

Nava, G., Fumagalli, F., Neutzner, S., and Di Fonzo, F. (2018). Large area porous 1d photonic
crystals comprising silicon hierarchical nanostructures grown by plasma-assisted, nanoparticle
jet deposition. Nanotechnology, 29(46):465603.

Ng, B. and Reid, W. (1979a). An initial value method for eigenvalue problems using compound

matrices. Journal of Computational Physics, 30(1):125–136.

Ng, B. and Reid, W. (1979b). A numerical method for linear two-point boundary-value problems

using compound matrices. Journal of Computational Physics, 33(1):70–85.

Ng, B. and Reid, W. (1985). The compound matrix method for ordinary differential systems.

Journal of Computational Physics, 58(2):209–228.

Nolte, A. J., Rubner, M. F., and Cohen, R. E. (2005). Determining the young’s modulus of polyelec-
trolyte multilayer films via stress-induced mechanical buckling instabilities. Macromolecules,
38(13):5367–5370.

Ogden, R. W. (1972). Large deformation isotropic elasticity–on the correlation of theory and
experiment for incompressible rubberlike solids. Proceedings of the Royal Society of London. A.
Mathematical and Physical Sciences, 326(1567):565–584.

Park, S., Wang, G., Cho, B., Kim, Y., Song, S., Ji, Y., Yoon, M. H., and Lee, T. (2012). Flexible

molecular-scale electronic devices. Nature Nanotechnology, 7(7):438–442.

Rivlin, R. (1948). Large elastic deformations of isotropic materials. i. fundamental concepts.
Philosophical Transactions of the Royal Society of London. Series A, Mathematical and Physical
Sciences, 240(822):459–490.

Rivlin, R. (1949). Large elastic deformations of isotropic materials. v. the problem of flexure.
Proceedings of the Royal Society of London. Series A. Mathematical and Physical Sciences,
195(1043):463–473.

Roccabianca, S., Bigoni, D., and Gei, M. (2011). Long wavelength bifurcations and multiple neutral
axes of elastic layered structures subject to finite bending. Journal of Mechanics of Materials
and Structures, 6(1):511–527.

72

Roccabianca, S., Gei, M., and Bigoni, D. (2010). Plane strain bifurcations of elastic layered struc-
tures subject to finite bending: Theory versus experiments. IMA Journal of Applied Mathematics
(Institute of Mathematics and Its Applications), 75(4):525–548.

Sadeghian, H., Yang, C. K., Goosen, J. F. L., Van Der Drift, E., Bossche, A., French, P. J.,
and Van Keulen, F. (2009). Characterizing size-dependent effective elastic modulus of silicon
nanocantilevers using electrostatic pull-in instability. Applied Physics Letters, 94(22):1–4.

Sigaeva, T., Mangan, R., Vergori, L., Destrade, M., and Sudak, L. (2018). Wrinkles and creases
in the bending, unbending and eversion of soft sectors. Proceedings of the Royal Society A:
Mathematical, Physical and Engineering Sciences, 474:2212.

Stafford, C. M., Harrison, C., Beers, K. L., Karim, A., Amis, E. J., VanLandingham, M. R., Kim,
H.-C., Volksen, W., Miller, R. D., and Simonyi, E. E. (2004). A buckling-based metrology for
measuring the elastic moduli of polymeric thin films. Nature materials, 3(8):545.

Stafford, C. M., Vogt, B. D., Harrison, C., Julthongpiput, D., and Huang, R. (2006). Elastic moduli

of ultrathin amorphous polymer films. Macromolecules, 39(15):5095–5099.

Steigmann, D. and Ogden, R. (1997). Plane deformations of elastic solids with intrinsic boundary
elasticity. Proceedings of the Royal Society of London. Series A: Mathematical, Physical and
Engineering Sciences, 453(1959):853–877.

Tahk, D., Lee, H. H., and Khang, D.-Y. (2009). Elastic moduli of organic electronic materials by

the buckling method. Macromolecules, 42(18):7079–7083.

Torres, J. M., Stafford, C. M., and Vogt, B. D. (2009). Elastic modulus of amorphous polymer
thin films: Relationship to the glass transition temperature. ACS Nano, 3(9):2677–2685. PMID:
19702280.

Triantafyllidis, N. (1980). Bifurcation phenomena in pure bending. Journal of the Mechanics and

Physics of Solids, 28(3-4):221–245.

Truesdell, C. and Noll, W. (2004). The non-linear field theories of mechanics. Springer.
Uchic, M. D., Dimiduk, D. M., Florando, J. N., and Nix, W. D. (2004). Sample dimensions

influence strength and crystal plasticity. Science, 305(5686):986–989.

Vendamme, R., Ohzono, T., Nakao, A., Shimomura, M., and Kunitake, T. (2007). Synthesis and
micromechanical properties of flexible, self-supporting polymer - sio2 nanofilms. Langmuir,
23(5):2792–2799. PMID: 17253729.

Wang, J., Zeng, Z., Weinberger, C. R., Zhang, Z., Zhu, T., and Mao, S. X. (2015). In situ atomic-
scale observation of twinning-dominated deformation in nanoscale body-centred cubic tungsten.
Nature Materials, 14(6):594–600.

Wang, L., Guan, P., Teng, J., Liu, P., Chen, D., Xie, W., Kong, D., Zhang, S., Zhu, T., Zhang, Z., Ma,
E., Chen, M., and Han, X. (2017a). New twinning route in face-centered cubic nanocrystalline
metals. Nature Communications, 8(1):1–7.

73

Wang, L., Teng, J., Sha, X., Zou, J., Zhang, Z., and Han, X. (2017b). Plastic Deformation through
Dislocation Saturation in Ultrasmall Pt Nanocrystals and Its in Situ Atomistic Mechanisms. Nano
Letters, 17(8):4733–4739.

Yue, Y., Liu, P., Zhang, Z., Han, X., and Ma, E. (2011). Approaching the theoretical elastic strain

limit in copper nanowires. Nano Letters, 11(8):3151–3155.

Zamiri, A. and De, S. (2011). Mechanical properties of hydroxyapatite single crystals from
nanoindentation data. Journal of the Mechanical Behavior of Biomedical Materials, 4(2):146–
152.

74