A tutorial on the stability and bifurcation analysis of the
electromechanical behaviour of soft materials
Shengyou Yang1,2
and Pradeep Sharma2,3∗
1
Department of Engineering Mechanics, School of Civil Engineering, Shandong University, Jinan, 250061,
China
2
Department of Mechanical Engineering, University of Houston, Houston, TX 77204, USA
3
Department of Physics, University of Houston, Houston, TX 77204, USA
∗
E-mail: psharma@uh.edu; Fax: +1-713-743-4503; Tel: +1-713-743-4502
Contents
1 Introduction 4
2 Stability and bifurcation: the big picture and basic concepts 7
2.1 The big picture . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 7
2.2 One-dimensional systems and geometric approach . . . . . . . . . . . . . . . . . . . . . . . 12
2.3 Deﬁnition of stability and the importance of norm . . . . . . . . . . . . . . . . . . . . . . 14
2.3.1 Stability analysis by linearization: one-dimensional examples . . . . . . . . . . . . 14
2.3.2 Principle of minimum energy: a conservative system . . . . . . . . . . . . . . . . . 15
2.3.3 Other stability criteria . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4 The concept of bifurcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.1 Bifurcation points . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 16
2.4.2 Bifurcation in a continuum media . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.5 More examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.5.1 Stability examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 18
2.5.2 Bifurcation examples . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 20
2.6 Closure on the basic concepts underlying stability and bifurcation . . . . . . . . . . . . . . 22
3 One-dimensional electrostatics of deformable dielectrics 22
3.1 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.1.1 Deformation of a body . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 22
3.1.2 Deformation gradient and strain . . . . . . . . . . . . . . . . . . . . . . . . . . . . 23
3.1.3 Diﬀerentiation and integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2 Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
1
Yang & Sharma: A tutorial on the stability and bifurcation Page 2
3.2.1 Spatial representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 24
3.2.2 Material representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 25
3.3 Balance of forces . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 26
3.4 Constitutive equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 27
3.5 Summary of the boundary-value problem in 1D . . . . . . . . . . . . . . . . . . . . . . . . 28
3.6 Incremental formulation and bifurcation analysis . . . . . . . . . . . . . . . . . . . . . . . 29
3.7 An energy formulation of the electrostatics of deformable dielectrics . . . . . . . . . . . . 30
3.7.1 Formulation in terms of the nominal electric ﬁeld ˜E1
. . . . . . . . . . . . . . . . . 30
3.7.2 Formulation in terms of the nominal polarization ˜P1
. . . . . . . . . . . . . . . . . 32
3.8 Examples of large deformation and bifurcation analysis in soft dielectrics . . . . . . . . . . 34
3.8.1 Large deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 35
3.8.2 Incremental boundary-value problem . . . . . . . . . . . . . . . . . . . . . . . . . . 35
4 Three-dimensional electrostatics of soft deformable dielectrics 37
4.1 Mathematical preliminaries . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.1.1 Basic tensor algebra . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 38
4.1.2 Tensor analysis: diﬀerentiation and integration . . . . . . . . . . . . . . . . . . . . 39
4.2 Kinematics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2.1 Deformation of a body . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 40
4.2.2 Deformation gradient, stretch and strain tensors . . . . . . . . . . . . . . . . . . . 40
4.2.3 Material and spatial gradient, divergence, and curl . . . . . . . . . . . . . . . . . . 42
4.2.4 Material and spatial integration . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 42
4.3 Maxwell’s equations . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3.1 Spatial representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 43
4.3.2 Material representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 44
4.4 Balance of forces and moments . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.4.1 Spatial representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 45
4.4.2 Material representation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 46
4.5 Constitutive equations of soft dielectrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.5.1 Decomposition of the stress tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.5.2 Elastic stress tensors . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 47
4.5.3 Maxwell stress tensor for an ideal dielectric elastomer . . . . . . . . . . . . . . . . 49
4.6 Constitutive law for soft dielectrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 50
4.7 Summary of the boundary-value problem . . . . . . . . . . . . . . . . . . . . . . . . . . . 51
4.8 Incremental formulation and bifurcation analysis . . . . . . . . . . . . . . . . . . . . . . . 51
4.9 Example 1: Pull-in instability of a dielectric elastomer ﬁlm by using the incremental method 53
4.9.1 Large deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 54
4.9.2 Incremental boundary-value problem . . . . . . . . . . . . . . . . . . . . . . . . . . 55
4.9.3 Solution of the incremental boundary-value problem . . . . . . . . . . . . . . . . . 57
2
Yang & Sharma: A tutorial on the stability and bifurcation Page 3
4.10 Example 2: Wrinkle surface instability of a dielectric elastomer by using the incremental
method . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.10.1 Large deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 58
4.10.2 Small deformation simpliﬁcation . . . . . . . . . . . . . . . . . . . . . . . . . . . . 60
4.10.3 Incremental boundary-value problem . . . . . . . . . . . . . . . . . . . . . . . . . . 61
4.10.4 Solution of the incremental boundary-value problem . . . . . . . . . . . . . . . . . 62
5 Electromechanical instability and bifurcation of soft dielectrics: pull-in instability 64
5.1 Pull-in instability of a dielectric elastomer ﬁlm . . . . . . . . . . . . . . . . . . . . . . . . 65
5.1.1 Stability analysis by using the Hessian approach . . . . . . . . . . . . . . . . . . . 65
5.1.2 Bifurcation analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 67
5.2 Avoiding the pull-in instability of a dielectric elastomer ﬁlm . . . . . . . . . . . . . . . . . 69
5.2.1 Stability analysis by using the energy method . . . . . . . . . . . . . . . . . . . . . 69
5.2.2 Bifurcation analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 71
6 Stability of the homogeneous deformation of soft dielectrics: an alternative energy
formulation 72
6.1 A dielectric ﬁlm . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 72
6.1.1 Compressible soft dielectrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.1.2 Incompressible soft dielectrics . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 74
6.1.3 Example: a ﬁlm of incompressible neo-Hookean dielectric . . . . . . . . . . . . . . 76
6.2 A dielectric disc . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.2.1 Large deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 78
6.2.2 Stability analysis and actuation strain . . . . . . . . . . . . . . . . . . . . . . . . . 82
6.2.3 Pull-in instability vs. electric breakdown . . . . . . . . . . . . . . . . . . . . . . . . 84
7 Electromechanical cavitation and the snap-through instability in soft hollow dielectric
spheres 85
7.1 Formulation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 85
7.1.1 Energy of the system . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 86
7.1.2 First variation of the energy functional . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.1.3 Spherical deformation . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 87
7.1.4 Ideal dielectric elastomer . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 89
7.1.5 Solution of the boundary-value problem . . . . . . . . . . . . . . . . . . . . . . . . 89
7.2 Electromechanical cavitation and the snap-through instability . . . . . . . . . . . . . . . . 90
7.2.1 Ball’s classic cavitation problem . . . . . . . . . . . . . . . . . . . . . . . . . . . . 90
7.2.2 Eﬀect of the electric ﬁeld on cavitation and snap-through instability . . . . . . . . 92
8 Post-buckling analysis: “Lyapunov-Schmidt-Koiter” approach 92
8.1 Bifurcation analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 93
8.2 Initial post-buckling analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 94
3
Yang & Sharma: A tutorial on the stability and bifurcation Page 4
8.3 Stability analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 96
8.4 Example: Buckling of Euler’s column subjected to an electric ﬁeld . . . . . . . . . . . . . 99
8.4.1 Linear bifurcation analysis . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
8.4.2 Energy stability criterion . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . . 101
8.4.3 Post-buckling analysis: LSK approach . . . . . . . . . . . . . . . . . . . . . . . . . 102
9 Concluding remarks and future study 104
Abstract
Soft materials, such as liquids, polymers, foams, gels, colloids, granular materials, and most
soft biological materials, play an important role in our daily lives. From a mechanical viewpoint,
soft materials can easily achieve large deformations due to their low elastic moduli; meanwhile,
surface instabilities, including wrinkles, creases, folds, and ridges, among others, are often observed.
In particular, soft dielectrics when subjected to electrical stimuli can achieve signiﬁcantly large
deformations that are often accompanied by instabilities. While instabilities are conventionally
thought to cause failures in the engineering context and carry a negative connotation, they can also
be harnessed for various applications such as surface patterning, giant actuation strain, and energy
harvesting. In the biological world, instability and bifurcation phenomena often precede important
events such as endocytosis, cell fusion, among others. Stability and bifurcation analysis (especially
for soft materials) is challenging and often presents a formidable barrier to entry in this important
ﬁeld. A multidisciplinary audience may lack the background in one or more areas that are needed to
carry out the requisite modeling or even understand papers in the literature. Furthermore, combining
electrostatics together with large deformations brings its own challenges. In this article, we provide
a tutorial on the basics of stability and bifurcation analysis in the context of soft electromechanical
materials. The aim of the article is to use simple examples and “gently” lead a reader, unfamiliar with
either stability analysis or electrostatics of deformable media, to develop the ability to understand the
pertinent literature that already exists and position them to embark on the state-of-the-art research
on this topic.
1 Introduction
The simple act of picking an object by a human-like robot would not be possible without the ability
of the underlying material to sustain large deformation. This example underscores the motivation
for studying soft materials. Beyond robotics [1–3] (Fig. 1), in contexts ranging from biotechnology to
electronics, soft materials are important wherever we need large deformations. Often, just the capacity
to sustain large deformation in response to mechanical forces is not enough. Ideally, we hope to
inﬂuence the material behavior by applying an external electric ﬁeld and this permits numerous applications
such as the already mentioned soft robotics, energy harvesting [4–7], sensors and actuators [8–11],
among others [12]. Soft dielectrics, when subjected to electrical stimuli, can achieve signiﬁcantly large
deformations (Fig. 2) due to the highly nonlinear electromechanical coupling between mechanical and
electric ﬁelds [13–15]. Large deformations invariably lead to the possibility of instabilities which include
the phenomena of surface wrinkles [16–22], creases [23–29], folds [30,31], electro-buckling [32–35], pull-in
4
Yang & Sharma: A tutorial on the stability and bifurcation Page 5
Figure 1: The utility of nonlinear behavior of soft materials. (a) A 9 cm tip-to-tip PneuNet actuator
gripping a raw chicken egg [1]. The actuator legs curled downward when the pressure is applied to the
PneuNets in the top active layer, and the soft gripper did not damage the egg. (b) A soft gripper made
of a buckling actuator that can pick up a piece of chalk [2].
instability [14,36–40], bursting drops in solid dielectrics [41], among others. Two types of electromechanical
instabilities observed in experiments are shown in Fig. 3: morphological instability of a single drop
in a solid dielectric polymer under electric ﬁelds and regular azimuthal wrinkles in a dielectric membrane.
While instabilities are generally often avoided in engineering applications as they are considered to be
a sign of, or precursor to, “failure”, as well reviewed recently [42,43], they may be exploited in a positive
manner. For example, Fig. 1(b) demonstrates the utility of the buckling of elastomer beams under pressure
to achieve a remarkably large actuation in soft machines [2]. In the neighborhood of the threshold
of instability, a small perturbation (in this case, pressure) can induce a large output—in this instance,
the buckling of pillars and the result that the air chambers collapse in Fig. 1(b). Other examples of the
utility of elastic instabilities include inducing negative Poisson’s ratio behavior in cellular solids [44], a
“smart” polymer adhesive based on surface wrinkling [45], the rapid closure of the Venus ﬂytrap through
the snap-buckling instability [46], the advanced functionality within shape changing materials [47]. More
examples of applications predicated on the buckling phenomenon can be found in the review [43].
A multidisciplinary audience may lack the background in one or more areas that are needed to carry
out the requisite modeling of instability and bifurcations in soft electromechanical solids—or even to
understand the literature in this ﬁeld. For example, most books that explain large nonlinear deformation
mechanics are diﬃcult for non-mechanicians to follow. Combining electrostatics, together with large
deformations, brings on its own challenges. Finally, while several good expositions [48–61] do exist on
stability and bifurcation analysis, some don’t necessarily take a tutorial approach to educate a multidisciplinary
audience while others may not cover all the pertinent topics that are central to the theme
5
Yang & Sharma: A tutorial on the stability and bifurcation Page 6
Figure 2: (a) Large deformation of a circular disc of VHB induced by an applied voltage on the upper
and bottom surfaces coated with electrodes of carbon grease [14]. (b) An elastomer gripper opens at a
voltage of 3.0 kV, which can grasp objects when the voltage is removed [8].
Figure 3: (a) Spherical drops of conductive liquid in solid dielectrics become unstable and burst under
suﬃciently high electric ﬁelds [41]. (b) Azimuthal wrinkles in dielectric elastomer sheets clamped on a
circular rigid frame [34].
of this paper. This tutorial article is intended to ﬁll this gap and will gently guide researchers, using
simple examples, to the basics of stability analysis and electrostatics of deformable media. Clearly, a
mere article (especially intended to be a tutorial) cannot hope to be either comprehensive or have the
rigor that leading experts in this ﬁeld expect. From that viewpoint, our aspirations are quite modest—we
expect the tutorial to simply be a convenient starting point to facilitate the reader’s ability to undertake
a study of the more specialized literature. In particular, since the focus is on elementary concepts, our
approach in this tutorial emphasizes analytically solvable illustrative examples (as opposed to discussion
of computational methods or problems that are only amenable computationally).
This paper is organized as follows. In Sec. 2, we discuss the basics concepts of stability and bifurcation,
along with some simple illustrative examples and their key pertinent features. In Sec. 3, we present onedimensional
(1D) nonlinear electroelasticity in continuum media. This section will be a useful starting
point for those who wish to decouple the conceptual understanding of electrostatics of deformable media
from the tensorial machinery necessary for three-dimensional (3D) continuum mechanics. In addition,
6
Yang & Sharma: A tutorial on the stability and bifurcation Page 7
the boundary-value problem and its incremental formulation are presented for the bifurcation analysis of
the electromechanical behaviour of soft materials. In Sec. 4, we present the 3D theory of electroelasticity
and formulate both the boundary-value problem and its incremental formulation. We ise pull-in stability
and electro-wrinkling as examples to illustrate the incremental method in electroelasticity. In Sec. 5, we
present the stability and bifurcation analysis of pull-in instability of soft dielectrics by using the energy
formulation presented by Zhao and Suo [38]. We use these examples to highlight the diﬀerence between
stability and bifurcation. In Sec. 6, we use an alternative energy formulation to study pull-in instability
of a dielectric ﬁlm and a disc. The literature is rife with “diﬀerent” formulations of electroelasticity.
While equivalent in some sense, these theories can often appear to be quite diﬀerent depending on the
independent variable used or the variational principle invoked. Liu [62,63] recently provided a detailed
overview of the various approaches to formulate electro-magneto-elasticity and taking a cue from him, we
have highlighted this aspect as well by presenting illustrative examples that use diﬀerent formulations. In
Sec. 7, we follow the energy formulation by Liu and present a study of electromechanical cavitation and
the snap-through instability of a hollow sphere in a dielectric elastomer. In Sec. 8, we brieﬂy introduce
the “Lyapunov-Schmidt-Koiter” (LSK) approach for the initial post-bifurcation analysis. We ﬁnally
conclude in Sec. 9 by highlighting topics not covered, or inadequately covered, in the tutorial, which the
reader may consider for further study.
2 Stability and bifurcation: the big picture and basic concepts
In this section, we introduce the key elements pertaining to stability and bifurcation analysis. This
subject is both broad and deep, and its nuances have been extensively covered from both a purely
mathematical viewpoint (its natural home) as well as in terms of discipline-speciﬁc applications that
range from biology, economics, physics, engineering, and others. We will, rather quickly, specialize all
our discussion primarily to nonlinear physical systems that are conservative (as opposed to dissipative).
This permits the use of powerful and well-established approaches like the calculus of variations to study
the stability. In the concluding section ( 9), where we discuss topics not covered in this article, we point
to the literature on stability and bifurcation analysis for the non-conservative case. The current section
is strongly inﬂuenced by the following books: [48–53].
2.1 The big picture
The vast majority of the diﬀerential equations that purport to describe nature (and often human made
problems) are nonlinear. In contrast, 99 % of the engineering curriculum (and most of the graduate curriculum)
is dedicated to studying linear systems. We would not be able to hear music without nonlinear
eﬀects in our ears; there will be no turbulence in ﬂuids and plasticity in solids without nonlinearities;
self-assembly of nano structures, phase transformations, superconductivity, laser device operation and
neuron interactions are but just a few examples of the manifestation of nonlinear phenomena. There
are good reasons why our curriculum is focused on linear systems. They are important to develop intuition
(—which can, however, slide straight out of the window once we discuss some nonlinear systems).
7
Yang & Sharma: A tutorial on the stability and bifurcation Page 8
Moreover, we can usually solve linear systems—even if not in closed-form, fairly easily using numerical
methods. However, what if we could understand the nature of the nonlinear solutions without actually
having to solve those nasty equations? What information can we extract from them? This is the subject
of stability analysis and is a fascinating ﬁeld of research in many ﬁelds including engineering and physics.
For example, some nonlinear systems will lead to bifurcations, i.e., the presence of multiple solutions. If
we are aware of the nature of the bifurcations, we can often extract interesting insights prior to the hard
work of numerically solving the equations and prevent ourselves from making mistakes in the computational
solutions of those problems. Before we go further into the technical content, it is worthwhile to
examine some of the statements we have just made in more depth.
Consider a simple example of a nonlinear system—motion of a pendulum [48]. It’s governing equation
is ¨y + g
L sin y = 0. This equation cannot be solved analytically, at least not in terms of elementary
functions. In undergraduate courses, the assumption of sin y ∼ y for |y(t)| 1 is usually made, thus
linearizing the equation. We can certainly be more sophisticated and use one of the perturbation methods
(to obtain an approximate analytical solution) and, therefore, perhaps even extend the range of the
solution to y(t) that is not too small. Is there an easier route, however, to extract the basic behavior
without actually solving the equations? This system is so simple that even our physical intuition appears
to be able to do better than the linear solution. For example, at low energy, the pendulum will simply
swing around some location. If high enough energy is imparted, the pendulum will simply whirl over the
top. It should be obvious that the linearized or perturbation solution cannot handle the “whirling over
the top” part of the problem.
Figure 4: Schematic of the nature of stability of three rigid balls under the inﬂuence of gravity: stable,
unstable, and neutral (from left to right).
Nonlinear equations (and hence by extension the physical phenomena they represent) often have multiple
solutions and so-called bifurcations. One question we may ask is when should we expect multiple
solutions and when are they “stable”? We have not yet deﬁned stability, but here is a simple example
that provide a nice intuitive idea—precise deﬁnition will be given later on. Consider rigid balls shown
in Fig. 4 under the inﬂuence of gravity. The ﬁrst ball (starting from left) is somewhat “stable”. By the
latter, we intuitively mean that under small vibrations or perturbations, the ball is likely to get a bit
8
Yang & Sharma: A tutorial on the stability and bifurcation Page 9
disturbed, but it will stay put in its position. This conforms to our intuition that the ball on the left
is in a stable state. In contrast, the middle ball may be considered to be unstable because any small
perturbation is likely to move it from its original position.
Figure 5: Schematic of the buckling of Euler’s column with pinned-pinned ends.
Let us elaborate further on multiple solutions and bifurcations. The latter refers to a qualitative
change in the stability behavior of a system as some parameter changes. For example, consider buckling
of a column in Fig. 5 under compression. As long as the “parameter” (in this case, axial force F) is
below a certain value Fcr, the column deforms axially, i.e., its length shortens. In order to ﬁnd the
shortening of this length, we have to solve the partial diﬀerential equations of nonlinear elasticity. Above
a certain value Fcr of the load, the equations of elasticity no longer admit a unique solution and allow for
multiple possibilities. As we know from our everyday experience, beyond Fcr, the column will “buckle”
and adopt either of the two equilibrium states—up or down. The critical load at which the beam buckles
is the bifurcation point. The key point is that multiple solutions are quite common in nonlinear systems;
however, stability and bifurcation analysis gives us information about how stable those solutions are and
under what conditions. Furthermore, bifurcation is related to how both the number and nature of the
stable solutions may change as some parameter is varied.
Here is an example from the textbook [50] about bifurcation and also underscores the kind of questions
we often ask about nonlinear systems
¨y + 2β ˙y − λy + y3
= 0, t > 0, (1)
where the superimposed ˙( ) = d()/dt denotes the diﬀerentiation with respect to time t, β and λ are
parameters. Equation (1) represents a Duﬃng type oscillator which has applications in biology and
electrical engineering. The steady-state or equilibrium1
equations are straightforward to obtain and are
the solutions of λy − y3
= 0. The value y = 0 is a solution regardless of the value of the parameter λ
and y = ±
√
λ for λ ≥ 0. We can graphically represent the so-called bifurcation diagram (see Fig. 6). In
1Steady-state conditions are sometimes referred to in the literature as equilibria. This is a somewhat ill-advised use of
the word but we think it is worthwhile to bring this subtlety to the reader’s attention.
9
Yang & Sharma: A tutorial on the stability and bifurcation Page 10
this diagram, (λ, y) = (0, 0) is the bifurcation point. In short, three distinct steady-state solutions are
possible. Which speciﬁc steady-solution will be realized in practice? Some insights can be obtained from
a direct numerical solution (see [50] for details).
Figure 6: Bifurcation diagram for steady-state solutions of Eq. (1).
An interesting example of bifurcation is given by Strogatz [48] regarding social dynamics. Assume
that a small segment of the population has a resolute belief in something, e.g., right of women to vote.2
Imagine four overlapping segments of the populations: (i) Group A: A group that holds resolutely to
the belief and nothing can change their mind, (ii) Group B is the group with the opposing views (iii)
This group agrees with A but its members are somewhat uncommitted. This means that if someone
argues against A then B becomes AB (a person who is ambivalent or sees both sides of the arguments,
(iv) this fourth group agrees with B but faced with suﬃcient arguments from A, will also join group
AB. Interestingly enough, a simple model of diﬀerential equations can be used to represent the dynamics
between the various groups. It can be shown that at a critical population of “true believers” (i.e. group
A) there is a bifurcation point and everyone will end up agreeing with A and may be used to explain
how revolutions or massive changes in society can be realized by a small (but fanatically committed)
sub-segment of a population ( 13 % in one speciﬁc model). Similar models also exist for population dynamics,
i.e., depending on the critical parameters (resources, number of mating pairs, etc.), a population
may shrink below some parameters or grow.
There are other examples that illustrate the importance of understanding nonlinearities and stability.
An aircraft is subjected to a variety of forces as it ﬂies. Beyond a critical speed and depending on how
the aircraft has been designed, the vehicle can end up in a sustained stable oscillation state. This can
be disastrous for the plan—as you can imagine. Indeed actual cases have been documented in the ﬁfties
and sixties where planes disintegrated in mid-air due to this phenomenon of “ﬂutter”. In materials science,
defect nucleations, the formation of microstructure, endocytosis, frequency entrainment in biology
2which, by the way, became national law in the USA in 1920 after much struggle and debate. Individual states began
to slowly ratify the 19th amendment over the years to follow. The last state to ratify this law was Mississippi—in 1984.
10
Yang & Sharma: A tutorial on the stability and bifurcation Page 11
and numerous other physical phenomena are consequences of nonlinear behavior and closely tied to the
concept of stability.
We now turn to a very important part of the “big picture”. Most books by mathematicians on nonlinear
systems tend to focus on the stability analysis of ﬁnite-dimensional systems. In other words, discrete
or ﬁnite degrees of freedom systems like the pendulum or an oscillator. These systems are governed
exclusively by a system of ordinary diﬀerential equations. On the other hand, engineers have usually
focused on inﬁnite-dimensional or continuous systems, e.g., theory of elasticity, Navier-Stokes equation
of ﬂuid mechanics, and so forth. The latter involves the study of partial diﬀerential equations. Why this
disconnect? There is a good reason for this. Most of the rigorous results in the stability analysis are
only available for ﬁnite degrees of freedom system, and accordingly, mathematicians have tended to focus
textbooks only in this setting. Research monographs and papers, of course, are a diﬀerent story. The
truth is that one must understand the stability notions dealt in the host of mathematics books (those
dealing with discrete systems) to develop some intuition and an understanding of results that can be
proven rigorously. Many of those methods and theorems, if not exactly, at least in spirit, carry over to
the continuous case also. Furthermore, a continuous system can often be reduced to a discrete system.
However, we want to caution the reader that the transfer of the results proven in most mathematics
books on ﬁnite degrees of freedom systems to continuous systems is not to be taken for granted. We
introduce a simple ﬁnite-dimensional system below along with some deﬁnitions.
Consider the spring-mass-dashpot system. This is essentially a damped oscillator (see Sec. 1.2 in the
textbook by Strogatz [48] for more details). The governing equation is
m¨y + b ˙y + ky = 0. (2)
This equation is linear and is said to be ﬁnite-dimensional with dimension 2. How do we come up
with this nomenclature? Note that we need two pieces of information to completely describe the system,
y and ˙y. This can be much easily seen if we express this equation in the standard format of nonlinear
systems, which involves expressing equations as coupled ﬁrst-order equations. To do this for Eq. (2), we
deﬁne two variables, y1 = y and y2 = ˙y. Then, we can convert Eq. (2) into
˙y1 = y2, ˙y2 = −
b
m
y2 −
k
m
y1. (3)
In general, the system can be n-dimensional and, of course, nonlinear:
˙y1 = f1(y1, y2, · · · , yn),
˙y2 = f2(y1, y2, · · · , yn),
...
˙yn = fn(y1, y2, · · · , yn).



(4)
11
Yang & Sharma: A tutorial on the stability and bifurcation Page 12
The key point is that in ﬁnite-dimensional systems, we are dealing with a countable ﬁnite number of
variables. We may express the above equations more concisely using vector and matrix notation:
˙y = f(y). (5)
Here, y(t) ∈ Rn
is an n-dimensional vector with components (y1, y2, · · · , yn)T
and f is a nonlinear
function of y. In the case of linear systems, we can reduce the system of equations to a matrix form:
˙y = My. (6)
Evidently, for the simple damped oscillator in Eq. (3), y(t) = (y1, y2)T
and M is


0 1
−k/m −b/m

.
In an inﬁnite-dimensional system, y is a ﬁeld and depends continuously in some n-dimensional regions.
Most inﬁnite-dimensional engineering problems are deﬁned in 3D regions (at most), so in that case, we
may say that y depends on both time and has spatial variation. In that case, f is an operator. Simply
put, in the case of continuous or inﬁnite-dimensional systems, we must contend with partial diﬀerential
equations such as the Navier-Stokes equation of ﬂuid mechanics or elasticity.
Strogatz’s book has a nice table (see Fig. 1.3.1 in [48]), which discusses various types of physical
systems in terms of dimensions and whether the systems are linear or not. The table also covers nonlinear
continuum theories of elasticity, plasticity, biological membranes, among others. We note that the
“coupled nonlinear oscillators” in the category of n 1 includes molecular dynamics.
2.2 One-dimensional systems and geometric approach
The 1D system is the simplest dynamical system and permits a rather simple introduction of several
essential concepts. We basically contend with systems like ˙y = f(y), where ˙y = dy/dt, y : I → R, I is an
open interval of R, and f : R → R. While in an actual physical problem, y(t) may not be position and t
may not be time, it is often useful to think of them in such terms as an aid to build intuition. In that
case, ˙y = dy/dt is the particle velocity whose motion is conﬁned to R. If we plot ˙y vs. y, then ˙y has the
character of a vector. At each point y, ˙y represents a vector with a direction.
Consider the following example from Strogatz [48]: ˙y = sin y. What are its equilibrium points or
ﬁxed points or steady-state solutions? It is simple to ﬁnd. We merely set ˙y = 0. This tells us that the
ﬁxed points of this system are y = nπ, where n ∈ Z are integers.
The plot of ˙y vs. y is called the phase portrait and represents how the trajectories ﬂow. Although
such plots may appear non-intuitive at ﬁrst, they provide a wealth of information without solving any
equations in both 1D and 2D; eventually establishing some level of intuition for higher dimensional sys-
tems.
12
Yang & Sharma: A tutorial on the stability and bifurcation Page 13
Figure 7: Phase portrait of the nonlinear diﬀerential equation ˙y = sin y. The arrows point to the right
when the velocity is positive ˙y(t) > 0 and to the left when ˙y(t) < 0.
Keeping in mind the spirit expressed in the preceding paragraph, we attempt to plot the phase diagram
without solving the nonlinear diﬀerential equation. This case, ˙y = sin y, is rather simple, of course,
and the result is shown in Fig. 7. The ﬁxed or equilibrium or steady-state solutions, y = nπ, are where
the curve intersects the y-axis.
Now we draw arrows to indicate the direction of the velocity vector. We point the arrow towards
the right when velocity is positive and vice-versa. We choose to draw the arrows on the y-axis in Fig. 7,
although it can be done anywhere. Notice that there are two types of ﬁxed points: those that have arrows
diverging from them and those which have arrows pointing towards them. Although we will deﬁne
stability more formally in due course, the ﬁgure graphically shows that stable points (marked by a solid
circle) are those towards which trajectories gravitate towards, and the reverse is true for unstable ﬁxed
points (marked by a hollow circle) in Fig. 7.
The phase portrait can be used to construct the actual solution qualitatively. To see this, consider
Fig. 7. Assume that our initial condition is y = 0. The phase portrait tells us that at that point, the
velocity is zero. In other words, the particle will not move. Let the initial condition be y = π/4. In the
beginning, the velocity y will increase due to ˙y > 0, then after reaching a maximum value, there will
be a period of deceleration ( ˙y > 0 but ˙y gradually decreases) until it approaches the ﬁxed point y = π.
Other trajectories, starting from diﬀerent initial points, can be similarly drawn. For more examples of
the phase portrait, for example, n-dimensional systems, the reader is referred to the textbooks [48–53].
13
Yang & Sharma: A tutorial on the stability and bifurcation Page 14
2.3 Deﬁnition of stability and the importance of norm
It is now time to deﬁne stability more formally. Although we have so for explicitly only given examples in
the context of 1D systems, the deﬁnition articulated here is valid for a general ﬁnite-dimensional system.
In general, the concept of stability is deﬁned by the motion of the physical system. Thus it is a dynamic
concept, and we have the following deﬁnition.
Deﬁnition of Stability: An equilibrium solution y0 at t = t0 of the system ˙y = f(y) is Lyapunov
stable if given a small number > 0, there exists δ( ) > 0 for any small perturbation, such that the
subsequent state y(t), t0 < t < ∞, satisﬁes y(t) − y0 < .
This implies that for any equilibrium point, any solution that is close to it initially (within some
regions δ of y0), remains so and is conﬁned within some regions of y0 for all positive time.
The concept of asymptotic stability requires that, eventually, any solution that starts close to the
equilibrium point must approach y0 as t → ∞ or y(t) − y0 → 0. Detailed mathematical discussion of
the deﬁnition of stability may be found in [64–68].
The issue of the norm and ﬁnite-dimensional systems vs. inﬁnite/continuous systems:
The symbol · is a suitable norm. This measures how close or far the two solutions are. The Euclidean
norm3
is frequently used in inﬁnite-dimensional systems; however, a very key result is that, for a ﬁnitedimensional
system, once a theorem has been proven in one norm, it also holds true in any other norm [69].
The proof of this very powerful statement is beyond the scope of the present article. Stated slightly
diﬀerently, we can say that in the context of a ﬁnite-dimensional system, if a solution is stable/unstable
in one particular norm, then it will be so in any other norm also. This result is not true for inﬁnitedimensional
or continuous systems. This, among other issues, is one of the key issues that make the
stability analysis of continuous systems more complicated. We refer the reader to Como and Grimaldi [69]
for a very nice example of how the change of norm can lead to diﬀering conclusions in a structural
mechanics problem.
2.3.1 Stability analysis by linearization: one-dimensional examples
Let us consider the case of a single diﬀerential equation, i.e., ˙y(t) = f(y, λ), where λ is the load parameter.
The Taylor expansion of f(y, λ) about a stationary state y = y0 at a ﬁxed load parameter λ reads
˙y(t) = f(y0, λ) +
∂f
∂y s
(y − y0) + o(|y − y0|). (7)
Here ∂f
∂y
s
denotes the partial derivative of f(y, λ) with respect to y at (y0, λ), i.e., ∂f
∂y
s
= ∂f
∂y (y0, λ).
The notation o(|y − y0|) denotes a function with the limiting behavior, i.e., the ratio |o|
|y−y0| approaches
3The Euclidean norm of a vector v = (v1, v2, ..., vn) ∈ Rn, for example, is deﬁned as v := v2
1 + v2
2 + ... + v2
n.
14
Yang & Sharma: A tutorial on the stability and bifurcation Page 15
zero as |y − y0| → 0. By dropping the zero term, f(y0, λ) = 0 for a stationary state y = y0, and
introducing a new variable h(t) = y(t) − y0, the linearized form of Eq. (7) is
˙h(t) =
∂f
∂y s
h(t). (8)
By inserting the ansatz h(t) = h0eµt
into Eq. (8), we can easily get µ = ∂f
∂y
s
. For a negative
µ = ∂f
∂y
s
< 0, the perturbed state h(t) = h0eµt
decreases as the time t increases, thus the stationary
state y = y0 is stable. Otherwise, the stationary state y = y0 is unstable for µ > 0 and a saddle point
for µ = 0.
Consider the nonlinear diﬀerential equation ˙y = sin y whose stability has already been investigated
by using the geometric approach in Sec. 2.2. The ﬁxed point is denoted by y0 = ±nπ, n ∈ Z, at
which sin y0 = 0. By Eqs. (7) and (8), the linearized form of ˙y = sin y at y = y0 is denoted by
˙h(t) = (cos y0) h(t), where h(t) = y(t) − y0. For example, at y0 = 0, ±2π, cos y0 = 1 > 0, and the
ﬁxed points y0 = 0, ±2π are unstable; in contrast, cos y0 = −1 < 0 at y0 = ±π and the ﬁxed points
y0 = ±π are stable. Their properties of stability are also shown in Fig. 7 by using the geometric approach.
Figure 8: Phase portraits of the diﬀerential equation ˙y(t) = y2
+ λ for diﬀerent load parameter λ. The
arrows in each ﬁgure point to the right when ˙y(t) > 0 and to the left when ˙y(t) < 0.
Consider another example; ˙y(t) = f(y, λ) = y2
+λ. The stationary state y = y0 satisﬁes f(y0, λ) = 0,
which implies three cases (see Fig. 8): (a) two stationary points y0 = ±
√
−λ for λ < 0; (b) one repeated
stationary point y0 = 0 for λ = 0; (c) no stationary point for λ > 0. By ∂f
∂y
s
in Eq. (8), i.e., ∂f
∂y
s
= 2y0,
we can conclude that: (a) the stationary point y0 = −
√
−λ is stable while y0 =
√
−λ is unstable for
λ < 0; (b) the stationary point y0 = 0 is a saddle point for λ = 0. The stationary points and their
stabilities are graphically depicted in Fig. 8. For further examples, see, the books [48,49].
2.3.2 Principle of minimum energy: a conservative system
The preceding sections highlight stability as a dynamical concept. A stability analysis is related to the
motion of the system under small perturbations, and the type of analysis illustrated so far is valid for a
general system, including both nonconservative and conservative ones.
15
Yang & Sharma: A tutorial on the stability and bifurcation Page 16
For a conservative system—which will become the primary focus of this tutorial, the stability analysis
of the motion of diﬀerential equations is identical to the discussion of the minimization of potential energy.
The proof can be found in many textbooks, for example, in the exposition by Triantafyllidis [65]. A state,
denoted by the general state variable x, is stable if its potential energy I(x) is lower than that of all the
neighborhood states, x + δx, such that
I(x) ≤ I(x + δx), (9)
where δx is suﬃciently small and · is a properly deﬁned norm for the space of deformation functions.
2.3.3 Other stability criteria
In contrast to the linear stability analysis and the principle of minimum energy, other stability criteria
are also often used [54, 55]. The adjacent equilibrium stability criterion, for example, asserts that a
primary equilibrium state becomes unstable when there are other nearby equilibrium states. A detailed
discussion of the diﬀerence between the energy method and the adjacent equilibrium stability criterion can
be found in a recent paper [70]. For nonconservative systems, the Lyapunov stability criterion [52,53,64]
can be used to discuss the stability of the motion of a dissipative system [55]. Furthermore, for the
initial post-buckling of continuum media, the “Lyapunov-Schmidt-Koiter” (LSK) approach is widely
used [55,61,65,67]. The LSK approach will be discussed later in this tutorial, and a simple example of
the initial post-buckling of a dielectric column will be given.
2.4 The concept of bifurcation
Let us assume that the state of a physical system can be represented by a generalized state variable u.4
The generalized load on the system is denoted by λ. The equation
f(u, λ) = 0 (10)
determines the stationary state of a physical system with the state variable u and the bifurcation parameter
λ.
2.4.1 Bifurcation points
We ﬁrst assume that, for a some initial range of the load parameter λ, there is only one solution (branch)
u0, the so-called principal branch, such that f(u0, λ) = 0 in Eq. (10). In Fig. 9, the principal branch
u0 bifurcates to a bifurcated branch u0 + v when the bifurcation parameter λ increases to λb that corresponds
to the state ub. The pair (ub, λb) is called a bifurcation point. For example, in the example of
the column buckling, the unbuckled state would be called the principal branch, while the buckled state
4This generalized state variable u could be either a scalar/tensor variable or a scalar/tensor-valued function, which
depends on the actual physical system. For the latter case, we have to use the functional derivative in the calculus of
variations [56,71–73].
16
Yang & Sharma: A tutorial on the stability and bifurcation Page 17
(a) (b)
Figure 9: Schematic of the bifurcation states of a physical system f(u, λ) = 0. u is the generalized state
variable of a physical system, and λ is the bifurcation (load) parameter. (a) A bifurcation diagram. (b)
A limit point.
is the bifurcated branch.
Now the questions are how to determine the existence of the bifurcated branch and to ﬁnd the bifurcation
point. To answer these questions, we have to introduce the implicit function theorem whose
discussion, especially in the context of interest, may be found in many books [49,52,56].5
Theorem (Implicit function theorem). Let S ⊆ R1+1
and write points in the set S as (u, λ) with
u ∈ R, λ ∈ R. Let f : R1+1
→ R be C1
on S. Let (u0, λ0) ∈ S with f(u0, λ0) = 0. If the derivative
∂f
∂u
(u0, λ0) = 0,
then there exists an r > 0 and a C1
function g : B(λ0; r) → R with g(λ0) = u0 such that
f(g(λ), λ) = 0 and
d g(λ)
d λ
= −
∂f
∂u
(u, λ)
−1 ∂f
∂λ
(u, λ)
u=g(λ)
for all λ ∈ B(λ0; r) = {λ ∈ R : |λ − λ0| < r}. 2
Here we only consider one equation f(u, λ) = 0 with two unknowns u and λ. The unknown u is called
the dependent variable while the unknown λ is called either the free variable or independent variable. We
express the unknown u in terms of λ, such that u = g(λ). The implicit function theorem gives us the
condition under which we can do this locally and yields a relation for computing d g(λ)
d λ .
In particular, if f(ub, λb) = 0 and ∂f
∂u (ub, λb) = 0, the point (ub, λb) is called the bifurcation point
5At this stage of the tutorial, we will illustrate the bifurcation theory by using the derivative of a function rather than
the functional derivative. This arrangement highly simpliﬁes the discussion and helps establish the basic ideas underlying
bifurcations.
17
Yang & Sharma: A tutorial on the stability and bifurcation Page 18
(Fig. 9) at which the uniqueness of the solution branch u = g(λ) does not hold.
2.4.2 Bifurcation in a continuum media
In this section, we brieﬂy discuss bifurcations in a continuum media. More classic bifurcation problems
of elasticity (e.g., buckling of a naturally straight rod, a whirling rod, a beam, a plate, a cylindrical shell,
a spherical cap, etc.) can be found in Chapter 5 of the book by Antman [56].
Let us consider a 1D continuum media, B1 = {X ∈ R : 0 ≤ X ≤ L}. The generalized state variable
of the 1D physical system is denoted by u, which is a C2
mapping u : B1 → Rn
, i.e., u ∈ C2
([0, L]; Rn
).
With a load parameter λ ∈ R, the stationary state of the physical system is represented by
f(u(X), λ) = 0. (11)
At a load parameter λ0, there exists a state u0, such that f(u0, λ0) = 0. Consider a small variation
δu(X) with suﬃciently small δu(X) in the neighborhood of u0(X). The expansion of the function
f(u, λ0) at (u0, λ0) reads
f(u0 + δu, λ0) = f(u0, λ0) +
∂f
∂u
(u0, λ0)δu + o(δu), (12)
where ∂f
∂u (u0, λ0) is called the (Fr´echet) derivative of f at (u0, λ0), and limδu→0
o(δu)
δu = 0. A much
weaker notation of a derivative is that of a directional derivative, i.e., the (Gˆateaux) derivative. In this
tutorial, we don’t distinguish the diﬀerence between these two derivatives. Thus, we take
∂f
∂u
(u0, λ0)δu ≡
d
dτ
f(u0 + τδu, λ0)
τ=0
= lim
τ→0
f(u0 + τδu, λ0) − f(u0, λ0)
τ
. (13)
By the implicit function theorem, a necessary condition for the bifurcation of f(u, λ) = 0 at (u0, λ0)
is that the (Fr´echet) derivative of f at (u0, λ0) is not invertible, namely
f(u0, λ0) = 0 and
d
dτ
f(u0 + τδu, λ0)
τ=0
= 0 for suﬃciently small δu . (14)
2.5 More examples
2.5.1 Stability examples
The following simple example illustrates the use of the principle of minimum energy to study the stability
of some conservative physical systems.
Example: Consider the stability of rigid balls in Fig. 4. We assume that the rigid ball can only move
on a smooth curve, y = y(x), 0 ≤ x ≤ a, see Fig. 10. Choose the horizontal axis as the reference at which
the rigid ball has zero potential energy. Then the potential energy of the ball subjected to the gravity
18
Yang & Sharma: A tutorial on the stability and bifurcation Page 19
Figure 10: The stability of a rigid ball at diﬀerent positions of the curve y = y(x) subjected to the
gravity. The potential energy of the ball at state x is given by I(x) = mgy(x). The 3D vision of the
rigid ball on a curved surface can refer to Fig. 4.
at state x is denoted by
I(x) = mgy(x), (15)
where m is the ball mass and g is the gravitational acceleration. By the stability condition (9) and the
potential energy (15), the stability of a rigid ball at state x implies that
y(x) ≤ y(x + δx), (16)
which, for arbitrary δx with suﬃciently small |δx|, gives
∂y
∂x
= 0 and
∂2
y
∂x2
≥ 0. (17)
Therefore, the left rigid ball is stable, the middle ball is unstable, and the right ball is neutrally
stable. In Fig. 10, we plot the three states and their stabilities corresponding to the minimization of the
potential energy (15). This stability example simply involves ﬁnding the extremum of a scalar function
with one scalar variable.
Note that the state variable x in this example is a scalar variable due to its ﬁnite-dimensional nature
(a single rigid ball). In contrast, continuous physical systems often involve scalar/tensor functions rather
than scalar/tensor variables as state variables. These functions, for example, include the displacement in
deformable solids, the electric displacement in electrostatics, and the temperature in thermodynamics,
etc. Thus, calculus of variations [56,71–73], is often used in the stability analysis of deformable solids.
19
Yang & Sharma: A tutorial on the stability and bifurcation Page 20
2.5.2 Bifurcation examples
We present a few examples in this sub-section for a better understanding of the concept of bifurcation—
especially the knowledge of nonuniqueness of solution and bifurcation points by using the implicit function
theorem. For more advanced examples on bifurcation analysis, the following literature may be
consulted: [48,52,56,57,60,74].
(a) (b)
Figure 11: Schematic of bifurcation diagrams. (a) Pitchfork type, from the equation λx − x3
= 0 in
Eq. (18). (b) From the equation λ2
+ x2
− 1 = 0 in Eq. (19).
Example: Consider the equation
f1(x, λ) = λx − x3
= 0, x, λ ∈ R. (18)
Apparently, (x, λ) = (0, R) is trivial solution, i.e., the principal branch. At point (x0, λ0) = (0, 0),
we have f1(0, 0) = 0 and ∂f1
∂x (0, 0) = 0, indicating that (x0, λ0) = (0, 0) is a bifurcation point at which
a unique solution branch (x, λ) = (0, R) does not hold. From the bifurcation diagram in Fig. 11(a), we
can ﬁnd that two bifurcated branches initiate at the bifurcation point (0, 0).
Example: Consider the equation
f2(x, λ) = λ2
+ x2
− 1 = 0, x, λ ∈ R. (19)
At two points (x, λ) = (0, −1) and (0, 1), both f2 = 0 and ∂f2
∂x = 0 are satisﬁed. The bifurcation
diagram of this example is plotted in Fig. 11(b).
Example: Consider Euler’s buckling in a continuum media (see Fig. 12(a)). The equilibrium equation
20
Yang & Sharma: A tutorial on the stability and bifurcation Page 21
(a) (b)
Figure 12: Buckling of Euler’s column with pinned-pinned ends. (a) Schematic of one model of buckling
and the coordinates used for analysis. u(X) is the angle between the tangent of the buckled column
and the horizontal axis. (b) Bifurcation diagram and the ﬁrst three modes of buckling. The straight
column corresponds to the principal branch that is denoted by the vertical axis, u = 0. The ﬁrst three
bifurcation points are (0, λ1), (0, λ2), and (0, λ3). Here the notation is λ = F
EI , and λ1 = π2
L2 , λ2 = 4π2
L2 ,
and λ3 = 9π2
L2 .
of the incompressible Euler’s column is [56,74]
f(u(X), λ) =
d2
u
dX2
(X) + λ sin u(X) = 0, 0 < X < L, (20)
and the boundary conditions are
du
dX
(0) =
du
dX
(L) = 0. (21)
The unbuckled column admits a principal solution u0(X) = 0 (see Fig. 12(a)). We now linearize the
above boundary-value problem (20) and (21). Consider a small variation δu = u . By using Eq. (14),
the resulting incremental forms of Eqs. (20) and (21) are
d2
u
dX2
+ λu = 0, 0 < X < L, (22)
and
du
dX
(0) =
du
dX
(L) = 0. (23)
Nonzero solutions of Eqs. (22) and (23) are
un(X) = An cos( λnX) at λn =
n2
π2
L2
, (24)
where n ∈ Z∗
is a non-zero integer and An is the amplitude. The bifurcation diagram and the ﬁrst three
modes of buckling are shown in Fig. 12(b).
21
Yang & Sharma: A tutorial on the stability and bifurcation Page 22
2.6 Closure on the basic concepts underlying stability and bifurcation
The issues of bifurcation in the equilibrium solutions of nonlinear solids and their stability are two important
but diﬀerent topics [57,75–78]. Generally speaking, stability theory mostly studies the response
of deformed bodies or structures under disturbance, while bifurcation theory mainly focuses on the nonuniqueness
of solutions in solid mechanics.
These two diﬀerent concepts and the underlying theories are linked. For example, the buckling of
Euler’s column corresponds to a bifurcation point at which the unbuckled deformation becomes unstable.
Other examples of the relation between the stability and the bifurcation include an inﬂated elastic cylinder
[77], a soap ﬁlm spanning a ﬂexible loop [78], and surface instability of a compressed half-space of
elastic materials [16,70], etc. In the examples mentioned above, the bifurcation points and the unstable
critical points coincide with each other. Such straightforward examples, however, are rather few, making
a complete understanding of the connection between the two issues diﬃcult.
Two landmark works in the mechanics community were carried out by Ericksen and Toupin [75] and
Hill [76] for the relation between stability and bifurcation. They showed that there was no bifurcation
from the trivial solution if the corresponding state of the trivial solution was stable, i.e., stability implies
uniqueness, but the converse does not follow.
3 One-dimensional electrostatics of deformable dielectrics
In this section, before proceeding to a 3D formulation, we present a 1D theory of electrostatics of
deformable media. The objective is to establish the key ideas of electroelasticity in a simple manner
without contending with the tedium of tensor algebra.
3.1 Kinematics
Kinematics in mechanics of deformable solids describes the motion of a body without considering the
forces, or any other ﬁelds (such as electrical) that cause the motion [79–81]. In this tutorial, we only
focus on static systems.
3.1.1 Deformation of a body
Consider a 1D body occupying the domain Ω1
R ⊆ R in 1D Euclidean space. The 1D body can be
identiﬁed with the occupied domain in some ﬁxed conﬁguration, which we term the reference conﬁguration.
The domain in a reference conﬁguration is denoted by Ω1
R ⊆ R, and the 1D body identiﬁed with
Ω1
R = {X ∈ R : Xa < X < Xb} is called the reference body. A point X ∈ Ω1
R is called a material point.
The boundary of Ω1
R is denoted by ∂Ω1
R = {X = Xa & Xb}.
22
Yang & Sharma: A tutorial on the stability and bifurcation Page 23
Consider a smooth function χ : Ω1
R → R that assigns to each material point X ∈ Ω1
R a point
x = χ(X) = X + u(X) ∈ Ω1
, (25)
where x is called the spatial point and u(X) is the displacement that relates the material point X ∈ Ω1
R
and the spatial point x ∈ Ω1
.
The function χ in Eq. (25) is invertible and there exists a unique function χ−1
: Ω1
→ R that assigns
to each spatial point x ∈ Ω1
a material point
X = χ−1
(x) ∈ Ω1
R. (26)
Mathematically, with the properties of inverse functions, each spatial point x ∈ Ω1
must correspond
to no more than one material point X ∈ ΩR and vice versa. Physically, a particle, denoted by X ∈ Ω1
R,
cannot break into two particles, and two diﬀerent particles, for example, X1, X2 ∈ Ω1
R, cannot occupy
the same spatial position during deformation.
3.1.2 Deformation gradient and strain
The (scalar) deformation gradient of 1D deformation is deﬁned by
F =
∂χ(X)
∂X
= 1 +
∂u(X)
∂X
> 0, (27)
where ∂u(X)/∂X is the displacement gradient. The inverse of the deformation gradient is
F−1
=
1
F
=
∂χ−1
(x)
∂x
. (28)
The inverse F−1
in Eq. (28) is related to the derivative of the inverse function that can be given by
the inverse function theorem, such that f−1
(x) = 1/f (X), where the prime denotes the derivative.
By using Eqs. (25), (27) and (28), the deformation gradient F and the inverse F−1
relate to the
reference length element dX and the current length element dx as
dx = FdX and dX = F−1
dx. (29)
In addition, we often use strain to describe the deformation. Consider the change in the squared
lengths dx2
− dX2
. Here dx2
(or the expression (dx)2
) is the squared distance of the two spatial points
while dX2
is the squared original distance of the two material points. We may introduce two scalars Es
and es, namely
dx2
− dX2
= 2Es dX2
= 2es dx2
, (30)
where Es is the so-called Green strain and es is the Euler strain. Substituting Eq. (29) into Eq. (30), we
23
Yang & Sharma: A tutorial on the stability and bifurcation Page 24
have
Es =
1
2
(F2
− 1) and es =
1
2
(1 − F−2
). (31)
Alternatively, the strains can be expressed in terms of the displacement gradient in Eq. (27):
Es =
∂u
∂X
+
1
2
∂u
∂X
2
and es =
1
2
(1 − (1 +
∂u
∂X
)−2
). (32)
In particular, for small deformation ∂u
∂X 1, if the high order terms in Eq. (32) are omitted, then
the Green and the Euler strains are equivalent, such that
Es = es =
∂u
∂X
. (33)
3.1.3 Diﬀerentiation and integration
The derivatives in the reference and current conﬁgurations are respectively called the material and spatial
derivatives. They are connected by the deformation gradient F in Eq. (27). The material and spatial
derivatives are
∂
∂X
and
∂
∂x
. (34)
Consider a smooth scalar ﬁeld, which is represented by φ(x) in the current conﬁguration (Ω1
) and
by φR(X) in the reference conﬁguration (Ω1
R). Using the chain-rule, we have the following relation
∂φR(X)
∂X
=
∂φ(x)
∂x
∂χ(X)
∂X x=χ(X)
= F
∂φ(x)
∂x x=χ(X)
or
∂φ
∂X
=
∂φ
∂x
∂χ
∂X
= F
∂φ
∂x
. (35)
The latter (35)2 is more concise, where we have dropped the arguments and the subscript.
In addition, with the length element relation dx = FdX in Eq. (29), we have the following relation
between the spatial and material integrals
Ω1
φ(x)dx =
Ω1
R
φR(X)F(X)dX or
Ω1
φdx =
Ω1
R
φFdX. (36)
3.2 Maxwell’s equations
3.2.1 Spatial representation
Maxwell’s equations of electrostatics govern the behavior of materials in the presence of charges or
when subjected to a potential diﬀerence. In 1D electrostatics, the Maxwell equations in the current
conﬁguration are (SI units convention)
E1
= −
∂ξ
∂x
and
∂D1
∂x
= ρ1
f , (37)
where E1
is the (true) electric ﬁeld, ξ is the electric potential, D1
is the (true) electric displacement
ﬁeld, and ρ1
f is the (free) charge density (free charge per unit length in the current conﬁguration). Free
24
Yang & Sharma: A tutorial on the stability and bifurcation Page 25
charges are those which exist and remain in the body regardless of any applied electrostatic boundary
conditions. We remark that Maxwell’s equations, as written in most physics and non-mechanics books,
are always in the current conﬁguration.
The electric boundary conditions (voltage-controlled) on x = xa and x = xb are given by
ξ(xa) = ξa and ξ(xb) = ξb. (38)
For physical reasons, we may introduce another ﬁeld, the polarization of a body and is deﬁned as:
D1
= ε0E1
+ P1
, (39)
where ε0 is the vacuum permittivity and P1
is the (true) polarization ﬁeld. The vacuum does not polarize
and therefore the polarization P1
is zero in the vacuum and the relation (39) reduces to D1
= ε0E1
.
The electric constitutive relation in the current conﬁguration that is most frequently used is the one
that linearly relates the electric displacement to the electric ﬁeld:
D1
= εE1
, (40)
where ε is the permittivity of the material. While the relation in vacuum is expected to be exact, the
above equation is a constitutive choice and is applicable to only linear dielectrics.
A direct consequence of Eqs. (39) and (40) for linear dielectrics is
P1
= (ε − ε0)E1
. (41)
3.2.2 Material representation
Since we have to couple elasticity with electrostatics, we often need to invoke Maxwell’s equations in the
reference conﬁguration also. They are (SI units convention)
˜E1
= −
∂ξ
∂X
and
∂ ˜D1
∂X
= ˜ρ1
f , (42)
where ˜E1
is the (nominal) electric ﬁeld, ξ is the electric potential, ˜D1
is the (nominal) electric displacement,
and ˜ρ1
f is the (free) charge density (free charge per unit length in the reference conﬁguration).
The electric boundary conditions (voltage-controlled) on X = Xa and X = Xb are given by
ξ(Xa) = ξa and ξ(Xb) = ξb. (43)
25
Yang & Sharma: A tutorial on the stability and bifurcation Page 26
Recalling the relation (35), the nominal electric ﬁeld ˜E1
is related to the true electric ﬁeld E1
by
˜E1
= −
∂ξ
∂X
= −
∂ξ
∂x
∂x
∂X
= E1
F. (44)
The nominal electric displacement ﬁeld in the reference conﬁguration is equal to the true electric
displacement ﬁeld, namely6
˜D1
= D1
, (45)
which leads to the relation
˜D1
= εF−1 ˜E1
. (46)
3.3 Balance of forces
The contact forces on the boundaries x = xa and x = xb in 1D formulation are
t1
(xa) = t1
xa
and t1
(xb) = t1
xb
. (47)
In 1D statics, balance of forces in the deformed body Ω1
= {x ∈ R : xa < x < xb} is
t1
xb
− t1
xa
+
xb
xa
b1
(x)dx = 0, (48)
where b1
(x) is the body force density (per unit length) in the current conﬁguration. The smooth function
t1
: Ω1
→ R can be interpreted as the force (or the stress) on the 1D body. We deﬁne that t1
(x) > 0
corresponds to an extended (spatial) point x while t1
(x) < 0 corresponds to a compressed (spatial) point
x. Using integration by parts, an alternative form of Eq. (48) is
xb
xa
∂t1
(x)
∂x
+ b1
(x) dx = 0. (49)
Since Eq. (49) holds for any Ω1
= {x ∈ R : xa < x < xb}, diﬀerential form or local form of balance
of forces is given by
∂t1
(x)
∂x
+ b1
(x) = 0. (50)
With the relation (36), the integral of the body force in Eq. (48) can be written in the reference
conﬁguration as
Xb
Xa
b1
R(X)F(X)dX =
xb
xa
b1
(x)dx or
Xb
Xa
b1
RFdX =
xb
xa
b1
dx. (51)
6Actually, the deﬁnition (45) comes from the conservation of the charge during deformation, and a detailed derivation
will be given later on when we present the general 3D formulation. Here we can show that xb
xa
∂D1
∂x
dx = xb
xa
ρ1dx =
Xb
Xa
˜ρ1dX =
Xb
Xa
∂ ˜D1
∂X
dX =
Xb
Xa
∂ ˜D1
∂x
∂x
∂X
dX. For linear dielectrics, a direct consequence of the deﬁnition (45) is ˜D1 =
D1 = εE1 = εF−1 ˜E1.
26
Yang & Sharma: A tutorial on the stability and bifurcation Page 27
With the traction boundary conditions
t1
R(Xa) = t1
Xa
and t1
R(Xb) = t1
Xb
, (52)
together with Eqs. (48) and (51), the global form of balance of forces in the reference conﬁguration
becomes
t1
Xb
− t1
Xa
+
Xb
Xa
b1
R(X)FdX = 0. (53)
By using integration by parts, similar to the procedure used in Eqs. (49) and (50), we have the local
form of balance of forces in the reference conﬁguration
∂t1
(X)
∂X
+ b1
0(X) = 0, (54)
where the subscript ‘R’ is dropped and b1
0(X) = b1
(X)F is the body force (density) in the reference
conﬁguration.
It is important to emphasize that the forces (or stresses) t1
R and t1
in the two conﬁgurations are the
same. In contrast, the body forces b1
0 and b1 are related by b1
0 = b1
F. In the way we have presented the
1D formulation, there is no change of the cross-sectional area, but the length elements are related by
dx = FdX.
3.4 Constitutive equations
In electroelasticity, the stress t1
can be decomposed into two parts:
t1
= te
+ tM
, (55)
where te
is the elastic stress and tM
is the so-called Maxwell stress (to be deﬁned shortly). Note that
the nominal and true stresses are the same in 1D continuum mechanics. The elastic stress is a function
of the deformation gradient:
te
= f1
(F). (56)
In one widely used constitutive choice (linear Hooke’s law), the elastic stress can be represented by
te
= C(F − 1) = C
∂u
∂X
, (57)
where C is the elastic modulus for small deformation.
The divergence of the Maxwell stress can be regarded as the electric body force. The readers can
refer to the works [82–84] on further discussion of other deﬁnitions of the Maxwell stress. The electric
body force comes from the concept of Lorentz force. For a charge distribution ρ1
f in an electric ﬁeld
E1
, the Lorentz force per unit volume is equal to ρ1
f E1
. By the identiﬁcation, ∂tM
/∂x = ρ1
f E1
, i.e., the
27
Yang & Sharma: A tutorial on the stability and bifurcation Page 28
divergence of the Maxwell stress can be related to the electric body force, the 1D Maxwell stress of linear
soft dielectrics in the current conﬁguration as7
tM
=
ε
2
(E1
)2
. (58)
It follows from Eqs. (37), (40) and (58) that ∂tM
/∂x = εE1
∂E1
/∂x = E1
∂D1
/∂x = ρ1
f E1
. A
more general derivation of the Maxwell stress in 3D from the balance of force is shown in Sec. 4.5.3.
By Eqs. (44) and (46), the 1D Maxwell stress (58) can also be written as either tM
= ε
2 (F−1 ˜E1
)2
or
tM
= 1
2ε ( ˜D1
)2
. Thus, the total stress t1
in Eq. (55) can be written as
t1
= f1
(F) +
ε
2
(F−1 ˜E1
)2
. (59)
3.5 Summary of the boundary-value problem in 1D
The boundary-value problem of the electrostatics of deformable dielectrics in the reference conﬁguration
can be compactly summarized as:
∂t1
(X)
∂X
+ b1
0(X) = 0
equilibrium equation
, ˜E1
= −
∂ξ
∂X
,
∂ ˜D1
∂X
= ˜ρ1
f
Maxwell’s equations
, ˜D1
= εF−1 ˜E1
relation (linear dielectrics)
,
t1
R(Xa) = t1
Xa
, t1
R(Xb) = t1
Xb
mechanical BCs
, ξ(Xa) = ξa, ξ(Xb) = ξb
electric BCs
,
t1
= f1
(F) +
ε
2
(F−1 ˜E1
)2
constitutive equation
, F = 1 +
∂u(X)
∂X
kinematics
.



(60)
In the above representation, we have chosen the constitutive law representative of a linear dielectric,
and the nominal electric ﬁeld is our state variable (in addition to the deformation). Other state variables
can also be used in the formulation, and, based on the review by Liu [63], we will elaborate on that later
in the context of 3D formulation. For the sake of brevity, we have only emphasized the force-traction
and voltage-controlled boundary conditions in Eq. (60).
In contrast to Eq. (60), the boundary-value problem in the current conﬁguration may be summarized
as:
∂t1
(x)
∂x
+ b1
(x) = 0
equilibrium equation
, E1
= −
∂ξ
∂x
,
∂D1
∂x
= ρ1
f
Maxwell’s equations
, D1
= εE1
relation (linear dielectrics)
,
t1
(xa) = t1
xa
, t1
(xb) = t1
xb
mechanical BCs
, ξ(xa) = ξa, ξ(xb) = ξb
electric BCs
,
t1
= f1
(es) +
ε
2
(E1
)2
constitutive equation
, es =
1
2
(1 − (1 +
∂u
∂X
)−2
) X=χ−1(x)
Euler strain
.



(61)
7While it is common to introduce the Maxwell stress in this manner (as evident from many works), we believe that its
emergence is most natural and clear by ﬁrst constructing a suitable energy functional and then using a variational principle.
The reader is referred to the works [62,63] of Liu where this is elaborated further.
28
Yang & Sharma: A tutorial on the stability and bifurcation Page 29
3.6 Incremental formulation and bifurcation analysis
Of interest here is the bifurcation analysis of the electromechanical behaviour of 1D deformable dielectrics.
We focus on the onset of bifurcation from a given trivial solution to other solutions of the
boundary-value problem in Sec. 3.5. Based on the implicit function theorem in Sec. 2.4, the equations of
equilibrium have a non-trivial solution bifurcating from the trivial solution only if their linearized forms
possess a non-zero solution. The linearized equations describe the response of 1D deformable dielectrics,
in a state of equilibrium, to inﬁnitesimal increments of the deformation and the electric ﬁeld (i.e., either
the electric ﬁeld or the electric displacement). We now linearize the above boundary-value problem (60).
The inﬁnitesimal increment of the deformation χ is denoted by u (X), i.e., χ(X) → χ(X) + u (X),
with suﬃciently small u (X) and ∂
∂X u (X) . Consider the variation of the deformation gradient
F = ∂χ
∂X . By the deﬁnition (13), the inﬁnitesimal increment of the deformation gradient is
F =
d
dτ
∂
∂X
(χ(X) + τu (X))
τ=0
=
∂u
∂X
. (62)
It is to be understood that the consequence of the change χ(X) → χ(X) + u (X) is F → F + F .
Similarly, the inﬁnitesimal increment of the inverse F−1
is
(F−1
) =
d
dτ
∂
∂X
(χ(X) + τu (X))
−1
τ=0
= −
∂χ(X)
∂X
−2
∂u
∂X
= −F−2
F . (63)
By the variation F → F + F with suﬃciently small F , the inﬁnitesimal increment (F−1
) can
be obtained as follows: (F−1
) = d
dτ (F + τF )
−1
τ=0
= −F−2
F . We can also construct the following
derivation. Consider the identity FF−1
= 1. Suppose that the deformation is changed from χ(X) to
χ(X)+τu (X). Diﬀerentiating both sides of the identity by τ at τ = 0, we obtain F F−1
+F(F−1
) = 0,
and, hence, that (F−1
) = −F−2
F .
The inﬁnitesimal increment of the electric potential ξ is denoted by ξ (X), i.e., ξ(X) → ξ(X)+ξ (X),
with suﬃciently small ξ (X) and ∂
∂X ξ (X) . The inﬁnitesimal increment of the nominal electric ﬁeld
˜E1
in Eq. (42) is expressed as
˜E =
d
dτ
−
∂
∂X
(ξ + τξ )
τ=0
= −
∂ξ
∂X
. (64)
29
Yang & Sharma: A tutorial on the stability and bifurcation Page 30
The inﬁnitesimal increment ˜D of the nominal electric displacement ˜D1
in Eq. (46) is8
˜D = (εF−1 ˜E1
) = ε(F−1
) ˜E1
+ εF−1 ˜E = εF−1
(−F−1 ˜E1
F + ˜E ). (65)
The inﬁnitesimal increment t of the total stress t1
in (59) is
t =
∂f1
(F)
∂F
F + (F−1 ˜E1
) ˜D . (66)
Now the incremental version of the boundary-value problem (60) can be represented by
∂t (X)
∂X
= 0,
∂ ˜D
∂X
= 0, t (Xa) = t (Xb) = 0, ξ (Xa) = ξ (Xb) = 0
incremental BCs
, (67)
along with Eqs. (62)−(66).
3.7 An energy formulation of the electrostatics of deformable dielectrics
As opposed to a “stress centric” approach described in the preceding section, as done in purely mechanical
continuum mechanics, we may also start by formulating free energy and then invoke a suitable variational
principle to arrive at the theoretical formulation. In the energy formulation, the deformation gradient F
can be chosen as one state variable related to the mechanical ﬁeld. The other state variable related to the
electric ﬁelds can be chosen to be either the nominal electric ﬁeld ˜E1
, the nominal electric displacement
˜D1
, or the nominal polarization ˜P1
. A rigorous formulation that presents a relative comparison of
the diﬀerent formulations may be found in Liu [63]. To proceed in the 1D formulation and simplify
the discussion, we will assume the cross-section of the 1D structure to be constant (unit area) that is
independent of the deformation.
3.7.1 Formulation in terms of the nominal electric ﬁeld ˜E1
Consider a smooth mapping χ : Ω1
R → R. We use (χ, ˜E1
) as the independent variables to formulate the
1D theory of electrostatics. By the identity ˜E1
= −∂ξ/∂X and the electric and mechanical boundary
conditions
ξ(Xa) = ξa, ξ(Xb) = ξb, t1
R(Xa) = t1
Xa
, t1
R(Xb) = t1
Xb
, (68)
we deﬁne the energy functional of the 1D system is deﬁned as
Fe[χ, ξ] =
Xb
Xa
W1
(F, ˜E1
) −
Xb
Xa
b1
0 · χ − (t1
R · χ)
Xb
Xa
, (69)
8The detailed derivation is as follows: ˜D = (εF−1 ˜E1) = d
dτ
ε ∂
∂X
(χ + τu )
−1
· − ∂
∂X
(ξ + τξ )
τ=0
=
ε d
dτ
∂
∂X
(χ + τu )
−1
· − ∂
∂X
(ξ + τξ )
τ=0
+ ε ∂
∂X
(χ + τu )
−1
· d
dτ
− ∂
∂X
(ξ + τξ )
τ=0
=
ε − ∂χ
∂X
−2
∂u
∂X
− ∂ξ
∂X
+ε ∂χ
∂X
−1
−∂ξ
∂X
= ε(F−1) ˜E1 +εF−1 ˜E . We may also arrive at the following derivation.
By the variations F−1 → F−1 + (F−1) and ˜E1 → ˜E1 + ˜E , we have ˜D = d
dτ
[ε(F−1 + τ(F−1) )( ˜E1 + τ ˜E )]
τ=0
=
ε(F−1) ˜E1 + εF−1 ˜E .
30
Yang & Sharma: A tutorial on the stability and bifurcation Page 31
where b1
0 is the body force and t1
R is the traction force. The equilibrium state is then determined by the
variational principle [63]
min
χ
max
ξ
Fe[χ, ξ]. (70)
• Variation of the electric ﬁeld
Let us suppose that ˜E1
is the function required to make the free energy Fe stationary. Consider the
variation ˜E1
→ ˜E1
+ τ ˜E1
1 . Here τ is a suﬃciently small parameter and ˜E1
1 is an arbitrary function that
satisﬁes all the kinematically admissible deformations. The variation of the electric potential ξ related to
the electric reads ξ → ξ + τξ1, where the arbitrary potential function ξ1 equals zero at both end-points,
i.e., ξ1(Xa) = ξ1(Xb) = 0.
The ﬁrst variation of the energy functional (69) with respect to the electric potential admits
∂Fe[χ, ξ + τξ1]
∂τ τ=0
=
Xb
Xa
∂W1
∂ ˜E1
∂ξ1
∂X
= ξ1
∂W1
∂ ˜E1
Xb
Xa
−
Xb
Xa
ξ1
∂
∂X
∂W1
∂ ˜E1
= 0. (71)
By ξ1(Xa) = 0 and ξ1(Xb) = 0 and the standard lemma of the calculus of variations, Eq. (71) gives
∂
∂X
∂W1
∂ ˜E1
= 0. (72)
• Variation of the deformation
Let us suppose that the deformation χ makes the free energy Fe stationary. Consider the variation
χ → χ + τχ1. Here χ1 is arbitrary and satisﬁes all the kinematically admissible deformations. The
electric quantities are assumed to be independent of the variation of the deformation. Then condition
for the vanishing of the ﬁrst variation of the energy functional (69) with respect to the deformation reads
∂Fe[χ + τχ1, ξ]
∂τ τ=0
=
Xb
Xa
∂W1
∂F
∂χ1
∂X
−
Xb
Xa
b1
0 · χ1 − (t1
R · χ1)
Xb
Xa
= 0, (73)
which, with integration by parts and the lemma of the calculus of variations, yields the equilibrium
equation
∂
∂X
∂W1
∂F
+ b1
0 = 0 (74)
and the natural (traction) boundary conditions
∂W1
∂F
= t1
R at X = Xa & Xb. (75)
We may further partition the energy function W1
(F, ˜E1
) of a dielectric elastomer as:
W1
(F, ˜E1
) = W1
elast(F) −
ε
2
F−1
( ˜E1
)2
. (76)
31
Yang & Sharma: A tutorial on the stability and bifurcation Page 32
It follows from Eqs. (72) and (76) that
∂
∂X
εF−1 ˜E1
= 0, (77)
which, by Eq. (46), is the Maxwell equation (42)2 in the absence of free charges. By Eq. (76), the total
stress is
∂W1
∂F
=
∂W1
elast
∂F
+
ε
2
(F−1 ˜E1
)2
= f1
(F) +
ε
2
(F−1 ˜E1
)2
, (78)
which is exactly the total stress shown in Eq. (59). A suitable mechanical constitutive law may be chosen
by specifying W1
elast.
3.7.2 Formulation in terms of the nominal polarization ˜P1
We once again follow Liu [63]. Consider a smooth mapping χ : Ω1
R → R. We use (χ, ˜P1
) as the
independent variables to formulate the 1D theory of electrostatics. By using the identity ˜E1
= −∂ξ/∂X,
the 1D Maxwell equation (42)2 in the reference conﬁguration in the absence of free charges is
∂
∂X
−ε0F−1 ∂ξ
∂X
+ F−1 ˜P1
= 0. (79)
The electric and mechanical boundary conditions are
ξ(Xa) = ξa, ξ(Xb) = ξb, t1
R(Xa) = t1
Xa
, t1
R(Xb) = t1
Xb
. (80)
The free energy of the electrostatic system is
Fp[χ, ˜P1
] = Fmech
[χ] + Felct
[χ, ˜P1
]. (81)
We claim that the equilibrium state is determined by the variational principle [63]
min
χ
min
˜P 1
Fp[χ, ˜P1
] (82)
subjected to the Maxwell equation (79) (which can be treated as a constraint) and the boundary conditions
(80). The mechanical part Fmech
[χ] in Eq. (81) can be written as
Fmech
[χ] =
Xb
Xa
W1
elast(F) −
Xb
Xa
b1
0 · χ − (t1
R · χ)
Xb
Xa
(83)
while the electric part Felct
[χ, ˜P1
] for linear dielectric elastomers is
Felct
[χ, ˜P1
] =
ε0
2
Xb
Xa
F−1 ∂ξ
∂X
2
+ ξ −ε0F−1 ∂ξ
∂X
+ F−1 ˜P1
Xb
Xa
+
Xb
Xa
F−1
( ˜P1
)2
2(ε − ε0)
. (84)
• Variation of the polarization Suppose that ˜P1
is the function required to make the free energy Fp
stationary. Consider the variation ˜P1
→ ˜P1
+ τ ˜P1
1 . Here τ is a suﬃciently small parameter and ˜P1
1 is
32
Yang & Sharma: A tutorial on the stability and bifurcation Page 33
an arbitrary function that satisﬁes all the kinematically admissible deformations. The variation of the
electric potential ξ related to the polarization reads ξ → ξ + τξ1, where the arbitrary potential function
ξ1 equals zero at both end-points. The variation of the 1D Maxwell equation (79) reads
∂
∂X
−ε0F−1 ∂ξ1
∂X
+ F−1 ˜P1
1 = 0. (85)
The ﬁrst variation of the free energy (81) with respect to the polarization is
∂Fp[χ, ˜P1
+ τ ˜P1
1 ]
∂τ
τ=0
= ε0
Xb
Xa
F−1 ∂ξ
∂X
∂ξ1
∂X
+
Xb
Xa
F−1 ˜P1 ˜P1
1
ε − ε0
+ ξ1 −ε0F−1 ∂ξ
∂X
+ F−1 ˜P1
+ ξ −ε0F−1 ∂ξ1
∂X
+ F−1 ˜P1
1
Xb
Xa
.
(86)
By ξ1(Xa) = ξ1(Xb) = 0 and Eq. (85), the third term on the RHS of Eq. (86) can be written as
ξ −ε0F−1 ∂ξ1
∂X
+ F−1 ˜P1
1
Xb
Xa
=
Xb
Xa
∂
∂X
ξ −ε0F−1 ∂ξ1
∂X
+ F−1 ˜P1
1
=
Xb
Xa
−ε0F−1 ∂ξ
∂X
∂ξ1
∂X
+ F−1 ∂ξ
∂X
˜P1
1 .
(87)
Thus, the vanishing of the ﬁrst variation (86) may be reformulated as
∂Fp[χ, ˜P1
+ τ ˜P1
1 ]
∂τ
τ=0
=
Xb
Xa
∂ξ
∂X
+
˜P1
ε − ε0
F−1 ˜P1
1 = 0. (88)
Since Eq. (88) must be satisﬁed for arbitrary ˜P1
1 , it is easy to see that we require
˜P1
= −(ε − ε0)
∂ξ
∂X
. (89)
A consequence of Eqs. (79) and (89) is
∂
∂X
−εF−1 ∂ξ
∂X
= 0, (90)
which is the same as Eq. (77). By substituting Eq. (89) into Eq. (84) and using integration by parts and
Eq. (90), we can show that
Felct
[χ, ˜P1
] = −
Xb
Xa
ε
2
F−1 ∂ξ
∂X
2
. (91)
It follows from Eqs. (83) and (91) that the free energy (81) can then be written as
Fp[χ, ˜P1
] =
Xb
Xa
W1
elast(F) −
ε
2
F−1 ∂ξ
∂X
2
−
Xb
Xa
b1
0 · χ − (t1
R · χ)
Xb
Xa
. (92)
• Variation of the deformation Suppose that the deformation χ makes the free energy Fp stationary.
Consider the variation χ → χ + τχ1. Here χ1 is arbitrary and satisﬁes all the kinematically admissible
deformations. The electric quantities are assumed to be independent of the variation of the deformation.
33
Yang & Sharma: A tutorial on the stability and bifurcation Page 34
Then the vanishing of the ﬁrst variation of the free energy (92) with respect to the deformation is
0 =
∂Fp[χ + τχ1, ˜P1
]
∂τ
τ=0
=
Xb
Xa
∂W1
elast
∂F
+ tM ∂χ1
∂X
−
Xb
Xa
b1
0 · χ1 − (t1
R · χ1)
Xb
Xa
=
∂W1
elast
∂F
+ tM
− t1
R · χ1
Xb
Xa
−
Xb
Xa
∂
∂X
∂W1
elast
∂F
+ tM
+ b1
0 · χ1,
(93)
where tM
is the Maxwell stress
tM
=
ε
2
F−1 ∂ξ
∂X
2
=
ε
2
(F−1 ˜E1
)2
(94)
that is the same as the one in Eq. (78) (or Eq. (59)). Since χ1 is an arbitrary function, the zero ﬁrst
variation (93) gives the equilibrium equation
∂
∂X
∂W1
elast
∂F
+ tM
+ b1
0 = 0 (95)
and the natural (traction) boundary conditions
∂W1
elast
∂F
+ tM
= t1
R at X = Xa & Xb. (96)
A third type of the energy formulation is also possible—in terms of the deformation χ and the nominal
electric displacement ˜D1
—and is omitted here but the reader can refer to the work by Liu [63].
3.8 Examples of large deformation and bifurcation analysis in soft dielectrics
Figure 13: Schematic of the deformation of a circular disc of soft dielectrics subjected to an applied
voltage Φ. Two layers of compliant electrodes are coated on the upper and bottom surfaces. The disc
deforms from thickness H to λ3H, and both the upper and bottom electrodes gain an electric charge of
magnitude Q.
34
Yang & Sharma: A tutorial on the stability and bifurcation Page 35
3.8.1 Large deformation
Consider a circular disc of a soft dielectric subjected to an applied voltage Φ (see Fig. 13). Recall the
BVP (60), the electric and mechanical boundary conditions are
ξ(0) = 0, ξ(H) = Φ, t1
R(0) = t1
R(H) = 0. (97)
We specialize to homogeneous deformations and the deformation gradient F = ∂χ
∂X = λ3H
H = λ3 is a
constant. The Maxwell equation in the BVP (60), together with the electric boundary conditions (97),
gives the solution of the potential ξ(X) = X
H Φ. Then the nominal electric ﬁeld is
˜E1
= −
∂ξ
∂X
= −
Φ
H
. (98)
It follows from Eqs. (98) that the Maxwell stress in the BVP (60) is
tM
=
ε
2λ2
3
Φ
H
2
. (99)
The total stress in (60) then is
t1
= f1
(λ3) +
ε
2λ2
3
Φ
H
2
, (100)
which is independent of the coordinate X due to a constant λ3. Thus, the mechanical boundary conditions
in Eq. (97) ﬁnally give the governing equation
f1
(λ3) +
ε
2λ2
3
Φ
H
2
= 0, (101)
which governs the deformation λ3 of the disc subjected to the applied voltage Φ.
3.8.2 Incremental boundary-value problem
In contrast to Eq. (97), the incremental electric and mechanical boundary conditions are
ξ (0) = ξ (H) = 0 and t (0) = t (H) = 0. (102)
The inﬁnitesimal increments of the potential and displacement are assumed to be ξ and u = λ3X,
where λ3 is suﬃciently small and independent of the coordinates. It follows from Eqs. (62) and (64) that
F = λ3 and ˜E = −
∂ξ
∂X
. (103)
Since F = λ3 and F = λ3 are independent of the coordinates, the increment ˜D in Eq. (65) and the
incremental Maxwell equation ∂ ˜D
∂X = 0 in Eq. (67) implies ∂ ˜E
∂X = 0. By Eq. (103), we have ∂2
ξ
∂X2 = 0,
which, together with the incremental electric boundary conditions (102)1, gives ξ (X) = 0. Thus,
˜E = 0 and ˜D = ελ−2
3
Φ
H
λ3, (104)
35
Yang & Sharma: A tutorial on the stability and bifurcation Page 36
where λ3 is the solution to the governing equation (101). The inﬁnitesimal increment t in (66) then is
t =
∂f1
(λ3)
∂λ3
λ3 − ελ−3
3
Φ
H
2
λ3. (105)
Thus, the incremental mechanical boundary conditions (102)2 give ∂f1
(λ3)
∂λ3
− ελ−3
3
Φ
H
2
λ3 = 0,
which requires
∂f1
(λ3)
∂λ3
− ελ−3
3
Φ
H
2
= 0, (106)
for nonzero λ3 and Eq. (106) is the necessary condition for bifurcation. We can also obtain Eq. (106)
directly from the governing equation (101) by using the implicit function theorem in Sec. 2.4.1. The
procedure is the same as that of the ﬁrst two bifurcation examples in Sec. 2.5.2. Let x = λ3 and
λ = ε
2
Φ
H
2
in Eq. (101) for notation simplicity, we have the reformulated governing equation f3(x, λ) =
f1
(x) + x−2
λ = 0. The condition of bifurcation is ∂f3(x,λ)
∂x = ∂f1
(x)
∂x − 2x−3
λ = 0, which is exactly the
same as Eq. (106). In the following, we will consider two kinds of constitutive laws.
• Linear constitutive law Consider the constitutive law (57). The governing equation (101) becomes
C(λ3 − 1) +
ε
2λ2
3
Φ
H
2
= 0 (107)
while the necessary condition (106) for bifurcation is
C − ελ−3
3
Φ
H
2
= 0. (108)
• Nonlinear constitutive law Assuming the neo-Hookean constitutive law to describe the mechanical
behavior of the solid, the strain-energy function W1
elast in Eqs. (95) and (96) can be written as W1
elast(λ3) =
µ
2 (λ2
3 +2λ−1
3 −3), where µ is the shear modulus for small deformation. By the elastic stress te
= f1
(λ3) =
∂W1
elast/∂λ3, the governing equation (101) becomes
µ(λ3 − λ−2
3 ) +
ε
2λ2
3
Φ
H
2
= 0 (109)
and the necessary condition (106) for bifurcation is
µ(1 + 2λ−3
3 ) − ελ−3
3
Φ
H
2
= 0. (110)
Regarding the governing equation (109), a similar equation can be found, see eq.(28) in the work [42]:
µ(λ3 − λ−2
3 ) +
ε
λ3
3
Φ
H
2
= 0. (111)
The diﬀerence between Eq. (109) and Eq. (111) is due to the simpliﬁed 1D model in our tutorial
which does not consider the deformation of the cross-section.
36
Yang & Sharma: A tutorial on the stability and bifurcation Page 37
Figure 14: Bifurcation of the equilibrium state,
√
x − x4 − λ = 0 in Eq. (112), of a ﬁlm of neo-Hookean
ideal dielectrics subjected to an applied electric ﬁeld. The load parameter λ here is the nominal electric
ﬁeld ˜E3 ε/µ in equation (29) by Zhao and Wang [42], and the state variable x is the stretch λ3 in the
thickness direction. The limit point, (1 − xb, λb) = (0.37, 0.69), is found by bifurcation analysis.
To make a direct comparison to the stability analysis in the work [42], we consider the equilibrium
equation (111), which is reformulated as
f4(x, λ) = x − x4 − λ = 0, (112)
where x = λ3 and λ = Φ
H ε/µ for notation simplicity. By using the implicit function theorem in
Sec. 2.4.1 or considering the statement below Eq. (106), the existence of the bifurcation point requires
∂f4
∂x
=
1
2
1 − 4x3
√
x − x4
= 0, (113)
giving x = (1
4 )1/3
= 0.63 that is exactly the critical point for the onset of pull-in instability in the
work [42]. The bifurcation diagram is plotted in Fig. 14.
4 Three-dimensional electrostatics of soft deformable dielectrics
In this section, we summarize nonlinear electroelasticity in the context of three-dimensional solids. Several
books may be consulted for foundations of mechanics [57, 79–81, 85–88]. For the nonlinear ﬁeld
theory of deformable dielectrics, the following articles are good references: [58,59,62,63,82,84,89–96].
37
Yang & Sharma: A tutorial on the stability and bifurcation Page 38
Table 1: A brief summary of properties of a second-order tensor S [79–81].
Name Symbol How obtained Notes
Components Sij = (S)ij Sij = ei · Sej for basis {ei} Sej = Sijei
Transpose ST
b · Sa = a · ST
b, ∀ a, b ∈ R3
ST
= Sjiei ⊗ ej
Trace tr S tr S = Sii (in sum) for {ei} tr S = tr ST
= S : I
Inverse S−1
SS−1
= S−1
S = I (ST
)−1
= (S−1
)T
Matrix
representation
[S] [S] =
















S11 S12 S13
S21 S22 S23
S31 S32 S33








Determinant det S or |S| det S = Su·(Sv×Sw)
u·(v×w) , for any basis
{u, v, w}
det S = det ST
,
det S−1
= (det S)−1
Principal
invariants
I1(S),
I2(S),
I3(S)
I1(S)=tr S, I2(S) = 1
2 [tr2
(S) −
tr(S2
)], I3(S) = det S
Cayley-Hamilton
equation
4.1 Mathematical preliminaries
4.1.1 Basic tensor algebra
A (second-order) tensor can be seen as a linear transformation from R3
space to R3
space. Thus a tensor
S is a linear mapping of vectors to vectors; that is, for a vector a ∈ R3
,
b = Sa ∈ R3
. (114)
The (Cartesian) components Sij of the tensor S are
Sij = ei · Sej, (115)
where the dot in Eq. (115) denotes the inner (dot) product between vectors. An alternative (component)
form of b = Sa in Eq. (114) is given by
bi = Sijaj. (116)
In this tutorial, the vector form (i.e., Eq. (114)) rather than the component form (i.e., Eq. (116))
of tensor algebra is often used. In the following, we just summarize some commonly used properties of
second-order tensors in Table 1.
38
Yang & Sharma: A tutorial on the stability and bifurcation Page 39
4.1.2 Tensor analysis: diﬀerentiation and integration
For a scalar function φ(x) and a vector function v(x) of vector x, their gradients are
grad φ(x) =
∂φ(x)
∂xi
ei and grad v(x) =
∂vi(x)
∂xj
ei ⊗ ej, (117)
respectively. In addition, the divergence of a vector ﬁeld v(x) and a tensor ﬁeld S(x) may be deﬁned by:
div v(x) =
∂vi(x)
∂xi
and div S(x) =
∂Sij(x)
∂xj
ei. (118)
Deformation is a mapping from material points X to spatial points x. Thus, it is necessary to
distinguish between the gradient (or the divergence) with respect to X and x. In this tutorial, we use
“ ” for the gradient and “Div” for the divergence with respect to material points X, for example,
φ(X) =
∂φ(X)
∂Xi
ei, v(X) =
∂vi(X)
∂Xj
ei ⊗ ej, Div v(X) =
∂vi(X)
∂Xi
. (119)
The expressions (117)−(119) may be simply regarded as the operations in the Cartesian coordinates.
Actually, if xi(or Xi) and {ei} are taken as the general coordinates and the corresponding basis, the
gradient and the divergence in the general coordinates also have similar forms as in Eqs. (117)−(119).
Speciﬁcally, vector operations in orthogonal curvilinear coordinates, such as cylindrical polar coordinates
and spherical polar coordinates, as well as other concepts in tensor analysis, can be found in textbooks
[85–88,97].
Consider a scalar function φ(S) of a tensor variable S. The derivative ∂φ(S)/∂S is the tensor function,
namely
∂φ(S)
∂S
=
∂φ(S)
∂Sij
ei ⊗ ej. (120)
The chain-rule reads
˙
φ(S) =
∂φ(S)
∂S
: ˙S, (121)
where S is a function of t and the dot denotes the derivative with respect to t. The chain-rule is a simple
and eﬃcient way of calculating the derivative ∂φ(S)/∂S.
One form of the divergence theorem reads
∂ΩR
SN dA =
ΩR
Div S dV, (122)
where ΩR ∈ R3
is a bounded region with smooth boundary ∂ΩR, S is a diﬀerentiable tensor ﬁeld over
ΩR, and N is the outward unit normal to ∂ΩR. An important identity from the divergence theorem
(122) is
∂ΩR
u · SN dA =
ΩR
(u · Div S + S : u) dV, (123)
where u is a diﬀerentiable vector ﬁeld over ΩR.
39
Yang & Sharma: A tutorial on the stability and bifurcation Page 40
4.2 Kinematics
4.2.1 Deformation of a body
A body may occupy a domain in Euclidean space. We may identify the body with the occupied domain
in some ﬁxed conﬁguration, which is called the reference conﬁguration. The domain in a reference conﬁguration
is denoted by ΩR, and the body identiﬁed with ΩR is called the reference body. A point X ∈ ΩR
is called a material point. The (suﬃciently regular) boundary of ΩR is denoted by ∂ΩR (see Fig. 15).
Consider a smooth function χ : ΩR → Ω that assigns to each material point X ∈ ΩR a point
x = χ(X) = u(X) + X ∈ Ω, (124)
where x is called the spatial point and u(X) is the displacement vector that relates the material point
X ∈ ΩR and the spatial point x ∈ Ω (see Fig. 15). The smooth function χ is the deformation and
Ω = χ(ΩR) is the domain occupied by the body at the current conﬁguration.
Figure 15: Reference and current conﬁgurations of a body.
4.2.2 Deformation gradient, stretch and strain tensors
The deformation gradient is deﬁned by
F = χ(X) = u + I, Fij =
∂xi
∂Xj
=
∂ui
∂Xj
+ δij, (125)
where δij is the Kronecker delta. The deformation gradient F is a key measure of the deformation and
other related quantities (the various forms of strain tensors) are derived from it.
The (volumetric) Jacobian, J(X), relates to the change in volume between the reference, dV , and the
40
Yang & Sharma: A tutorial on the stability and bifurcation Page 41
current volume element, dv, as
dv = J(X)dV, (126)
where the Jacobian J(X) is equal to the determinant of the deformation gradient:
J(X) = det F(X) = det χ(X) > 0. (127)
In particular, the Jacobian for incompressible materials equals unity, i.e.,
J(X) = 1. (128)
The polar decomposition reads
F = RU = VR, (129)
where R is a rotation, U and V are positive-deﬁnite symmetric tensors. The explicit representations of
U and V are
U =
√
FT F and V =
√
FFT . (130)
Note that U and V are the so-called right and left stretch tensors, namely
λ = |U(X)e|, λ2
= e · U2
(X)e = e · C(X)e, (131)
where λ is the stretch at material point X in the material direction e, and C = U2
= FT
F is the right
Cauchy-Green tensor. In contrast, B = V2
= FFT
is the left Cauchy-Green tensor.
It is known that the eigenvalues of a positive-deﬁnite symmetric matrix (tensor) are real and positive,
and the eigenvectors corresponding to diﬀerent eigenvalues are orthogonal [97, 98]. Therefore, spectral
representations of the symmetric and positive-deﬁnite tensors U, V, C, and B are [79,80]
U =
3
i=1
λiri ⊗ ri, V =
3
i=1
λili ⊗ li, C =
3
i=1
λ2
i ri ⊗ ri, and B =
3
i=1
λ2
i li ⊗ li, (132)
where λi > 0 are the principle stretches, ri and li are the right and left principal directions.
In addition to the stretch, the strain tensors are often used to describe the deformation. Consider the
change in the squared lengths |dx|2
− |dX|2
. Here |dx|2
= dx · dx is the squared distance of two spatial
points in the current conﬁguration while |dX|2
is the squared original distance of two material points in
the reference conﬁguration. We introduce two second-order tensors Es
and es
, namely
|dx|2
− |dX|2
= 2dX · Es
dX = 2dx · es
dx, (133)
41
Yang & Sharma: A tutorial on the stability and bifurcation Page 42
where Es
is called the Green-St. Venant (Lagrangian) strain tensor or simply the Green strain tensor
Es
=
1
2
(FT
F − I) =
1
2
(C − I) =
1
2
( u + uT
+ uT
u) (134)
and es
is called the Almansi-Hamel (Eulerian) strain tensor or simply the Euler strain tensor
es
=
1
2
(I − F−T
F−1
) =
1
2
(I − B−1
). (135)
In linear elasticity, for small deformation u 1, the Green strain tensor Es
and the Euler strain
tensor es
are equivalent.
4.2.3 Material and spatial gradient, divergence, and curl
In the reference conﬁguration, the (material) gradient, the (material) divergence, and the (material) curl
are denoted, respectively, by
, Div, and Curl, (136)
while in the current conﬁguration, the (spatial) gradient, the (spatial) divergence, and the (spatial) curl
are denoted, respectively, by
grad, div, and curl. (137)
For a scalar ﬁeld ϕ and a vector ﬁeld g, with the chain-rule, we have the following relations between
the material and the spatial gradient:
∂ϕ
∂Xi
=
∂ϕ
∂xj
∂χj
∂Xi
= Fji
∂ϕ
∂xj
and
∂gi
∂Xj
=
∂gi
∂xk
∂χk
∂Xj
=
∂gi
∂xk
Fkj (138)
or, the vector forms
ϕ = FT
grad ϕ and g = (grad g)F. (139)
An equivalent form of Eq. (139) is
grad ϕ = F−T
ϕ and grad g = ( g)F−1
. (140)
In addition, the relation between the material and the spatial divergence reads:
Div g = FT
: grad g and div g = F−T
: g. (141)
4.2.4 Material and spatial integration
Recall Eq. (126), the spacial volume element dv is related to the material volume element dV through
the Jacobian J = det F, that is, dv = J dV . Consider a scalar ﬁeld which is represented by φ(x) in the
current conﬁguration (Ω) and by φR(X) in the reference conﬁguration (ΩR). With the volume element
42
Yang & Sharma: A tutorial on the stability and bifurcation Page 43
relation (126), dv = J dV , we have the following volume integral
Ω
φ(x)dv(x) =
ΩR
φR(X)J(X)dV (X) or
Ω
φdv =
ΩR
φJdV. (142)
The latter equation (142)2 admits a more concise form in which the arguments are suppressed and the
subscript of φ is dropped.
In contrast to the volumetric Jacobian J = detF in Eq. (126), the areal Jacobian Ja is deﬁned by
Ja = J|F−T
N|. (143)
Thus the spatial area element da = |n da| and the material area element dA = |N dA| have the relation
da = JadA, n da = JF−T
N dA = FC
N dA, (144)
where FC
= JF−T
is called the cofactor of F.
Consider a vector ﬁeld u and a tensor ﬁeld G. With Eq. (144), we have the following surface integrals
∂Ω
u · nda =
∂ΩR
u · FC
N dA and
∂Ω
Gnda =
∂ΩR
GFC
N dA. (145)
4.3 Maxwell’s equations
4.3.1 Spatial representation
In the absence of currents and magnetic ﬁelds, the Maxwell equations in the current conﬁguration become
(SI units convention)
curl E = 0 and div D = ρf , (146)
where E is the (true) electric ﬁeld, D is the (true) electric displacement, and ρf is the (free) charge
density in the current conﬁguration.
Equation (146)1 implies that the electric ﬁeld E in electrostatics is irrotational [97], that is, there
exists a scalar potential ﬁeld ξ, such that
E = −grad ξ. (147)
Equation (146)2 means that the electric ﬂux leaving a domain is proportional to the charge inside.
The electric displacement is deﬁned as
D = ε0E + P, (148)
43
Yang & Sharma: A tutorial on the stability and bifurcation Page 44
where P is the (true) polarization ﬁeld in the current conﬁguration. The relation (148) becomes D = ε0E
in vacuum. The constitutive relation for dielectrics can be written as D = ε(x)E, where ε(x) is the
second-order dielectric permittivity tensor. We will mostly study isotropic systems so ε(x) = εI, where
ε is a constant and I is the second-order identity tensor. The constitutive relation for linear dielectrics
in the current conﬁguration is simply:
D = εE. (149)
A direct consequence of Eqs. (148) and (149) for linear dielectrics is
P = (ε − ε0)E. (150)
4.3.2 Material representation
In contrast to Eq. (146), the Maxwell equations in the reference conﬁguration are
Curl ˜E = 0 and Div ˜D = ˜ρf , (151)
where ˜E is the nominal electric ﬁeld, ˜D is the nominal electric displacement, and ˜ρf is the (free) charge
density in the reference conﬁguration.
The irrotational ﬁeld ˜E in Eq. (151)1 indicates a scalar potential ﬁeld ξ, together with Eqs. (139)
and (147), we have
˜E = − ξ = −FT
grad ξ = FT
E. (152)
Recall the surface integral (145) of the spatial and material representations. If we replace the vector
ﬁeld u in Eq. (145) by the (true) electric displacement D, we have
∂Ω
D · nda =
∂ΩR
JD · F−T
N dA =
∂ΩR
JF−1
D · N dA =
∂ΩR
˜D · N dA, (153)
which indicates the relation between the nominal and true electric displacements, namely
˜D = JF−1
D. (154)
The 3D relation (154) can be directly reduced to the 1D relation (45) since the Jacobian J and deformation
gradient F in 1D formulation admit J = F.
The deﬁnition of the relation between the nominal and true polarizations is not unique. For example,
the deﬁnition in the work [63] is
˜P = JP. (155)
By Eqs. (148), (154) and (155), the auxiliary ﬁeld in the reference conﬁguration can be written as
˜D = ε0JC−1 ˜E + F−1 ˜P. (156)
44
Yang & Sharma: A tutorial on the stability and bifurcation Page 45
Table 2: Conversion formulas
Nominal/true electric displacement ( ˜D/D), electric ﬁeld (˜E/E), polarization (˜P/P)
˜D = ˜E = ˜P = D = E = P =
( ˜D, F) ˜D
C ˜D
εJ
ε − ε0
ε
F ˜D
F ˜D
J
F ˜D
εJ
ε − ε0
εJ
F ˜D
(˜E, F) εJC−1 ˜E ˜E (ε − ε0)JF−T ˜E εF−T ˜E F−T ˜E (ε − ε0)F−T ˜E
(˜P, F)
εF−1 ˜P
ε − ε0
FT ˜P
(ε − ε0)J
˜P
ε˜P
(ε − ε0)J
˜P
(ε − ε0)J
˜P
J
(D, F) JF−1
D
FT
D
ε
(ε − ε0)J
ε
D D
D
ε
ε − ε0
ε
D
(E, F) εJF−1
E FT
E (ε − ε0)JE εE E (ε − ε0)E
(P, F)
εJF−1
P
ε − ε0
FT
P
ε − ε0
JP
εP
ε − ε0
P
ε − ε0
P
The constitutive relation of linear dielectrics in the reference conﬁguration, for example, is
˜D = εJC−1 ˜E. (157)
Combining Eqs. (156) and (157) for linear dielectrics, we have the relation
˜P = (ε − ε0)JF−T ˜E. (158)
For convenience, several important relations related to electrostatics of a linear isotropic dielectric
deformable media are summarized in Table 2.
4.4 Balance of forces and moments
The mechanical laws, including balance of forces and moments, can be represented respectively in the
current and reference conﬁgurations. We ﬁrst discuss the balance laws and the stress tensor ﬁeld in the
current conﬁguration.
4.4.1 Spatial representation
The balance of forces and moments for the static equilibrium of a deformed body Ω are
∂Ω
t(n) da +
Ω
b0 dv = 0, (global) balance of forces (159)
45
Yang & Sharma: A tutorial on the stability and bifurcation Page 46
and
∂Ω
r × t(n) da +
Ω
r × b0 dv = 0, (global) balance of moments, (160)
where t(n, x) is the (true) traction force on the boundary ∂Ω, n(x) is the outward unit normal to the
boundary ∂Ω, b0(x) is the (true) body force in the domain Ω, and r(x) is the position vector.
A well-known theorem in continuum mechanics is Cauchy’s theorem for the existence of stress, that
is, there exists a spatial tensor ﬁeld σ, known as the Cauchy stress (ﬁeld), due to the balance of forces
Eq. (159), namely
t(n) = σn, ti = σijnj. (161)
In contrast to Eqs. (159) and (160), we have local forms of the force and moment balances by using
Cauchy’s theorem (161), the divergence theorem (122), the identity (123), and the property of skew
tensors, such that
div σ + b0 = 0, σij,j + b0i = 0, (local) balance of forces (162)
and
σ = σT
, σij = σji, (local) balance of moments. (163)
Clearly, the local form (163) is that the Cauchy stress tensor σ is symmetric. It should be noted that
Eqs. (161)−(163) are consequences of Eqs. (159) and (160).
4.4.2 Material representation
Recall the surface integral (145) of the spatial and material representations. If we replace the tensor ﬁeld
G in Eq. (145) by the Cauchy stress tensor σ, we have
∂Ω
σnda =
∂ΩR
σFC
N dA =
∂ΩR
TN dA, (164)
where
T = σFC
= JσF−T
(165)
representing the stress measured per unit area in the reference conﬁguration. In contrast to the Cauchy
stress σ in Eq. (161), T in Eq. (165) is referred to as the (ﬁrst) Piola-Kirchhoﬀ stress.
With the volume integral (142), we obtain
Ω
b0 dv =
ΩR
J b0 dV =
ΩR
b0R
dV, (166)
where
b0R
= J b0, (167)
46
Yang & Sharma: A tutorial on the stability and bifurcation Page 47
is the body force measured per unit volume in the reference conﬁguration.
Combining Eqs. (164) and (166) and Cauchy’s theorem (161), the spatial form of balance of forces
(159) is equivalent to
∂ΩR
TN dA +
ΩR
b0R
dV = 0, (global) balance of forces. (168)
With the divergence theorem (122) in an arbitrary chosen domain ΩR, we have the local form
Div T + b0R
= 0, Tij,j + b0Ri = 0, (local) balance of forces. (169)
An alternative form of the stress relation (165) is
σ = J−1
TFT
. (170)
Since the Cauchy stress is symmetric, σ = σT
in Eq. (163), we have the material form of balance of
moments
TFT
= FTT
, TikFjk = FikTjk, (local) balance of moments. (171)
4.5 Constitutive equations of soft dielectrics
4.5.1 Decomposition of the stress tensors
In the electrostatic of soft dielectrics, both the (ﬁrst) Piola-Kirchhoﬀ stress T in Eq. (165) and the
Cauchy stress σ in Eq. (170) can be decomposed into two parts:
T = Te
+ TM
(172)
and
σ = σe
+ σM
. (173)
Here Te
and σe
are elastic stresses. In contrast, TM
in Eq. (172) is the Piola-Maxwell stress and σM
in Eq. (173) is the true Maxwell stress
4.5.2 Elastic stress tensors
Consider a strain-energy function We
= We
(F). Evidently, We
= We
(F) is a scalar-valued function of
one tensor variable F. The elastic part Te
of the Piola-Kirchhoﬀ stress in Eq. (172) is
Te
=



∂We
(F)
∂F
∂We
(F)
∂F
− LaF−T
, Te
ij =



∂We
(Fij)
∂Fij
, compressible solids,
∂We
(Fij)
∂Fij
− LaF−1
ji , incompressible solids,
(174)
47
Yang & Sharma: A tutorial on the stability and bifurcation Page 48
where La is the Lagrange multiplier related to the constraint of incompressibility J = det F = 1. Also,
La in Eq. (174) may be interpreted as the hydrostatic pressure. With Eqs. (170) and (174), the elastic
part σe
of the Cauchy stress in Eq. (173) is
σe
= J−1
Te
FT
, σe
ij = J−1
Te
ikFjk. (175)
In the following, we will brieﬂy discuss the consequences of frame-indiﬀerence in the strain-energy
functions and stresses. Consider a rotation Q ∈ Orth+
= {all rotations}. The frame-indiﬀerence of the
strain-energy functions reads [80]
We
(F) = We
(QF). (176)
With the polar decomposition F = RU, R ∈ Orth+
, in Eq. (129), and letting Q = RT
in Eq. (176), we
have
We
(F) = We
(RT
RU) = We
(U). (177)
Since U =
√
UT U =
√
UT RT RU =
√
FT F =
√
C, we can introduce a strain-energy function
¯We
(C), such that
¯We
(C) = We
(
√
C) = We
(U) = We
(F). (178)
Similarly, with the Green strain tensor Es
= 1
2 (C − I), Eqs. (177) and (178), we can have a strainenergy
function ˜We
(Ee
), namely
˜We
(Ee
) = ¯We
(C) = We
(U) = We
(F). (179)
The derivatives of the strain-energy functions in Eq. (179) have the following relations
∂We
(F)
∂F
= 2F
∂ ¯We
(C)
∂C
= F
∂ ˜We
(Ee
)
∂Ee
. (180)
In addition, assuming isotropy, the strain-energy function We
depends on F through the principal
stretches λ1, λ2 and λ3:
We
(F) = We
(λ1, λ2, λ3). (181)
In contrast to Eqs. (174) and (175), the elastic part of the ﬁrst Piola-Kirchhoﬀ stress and the Cauchy
stress can be written as
Te
=
3
i=1
∂We
∂λi
−
1
λi
La ˆni ⊗ ˆNi (182)
and
σe
=
3
i=1
J−1
λi
∂We
∂λi
− La ˆni ⊗ ˆni, (183)
where the orthonormal vectors ˆni and ˆNi are the spatial and referential directions, respectively; and
J = λ1λ2λ3 and La is zero for compressible solids while J = 1 and La is the hydrostatic pressure for
incompressible solids.
48
Yang & Sharma: A tutorial on the stability and bifurcation Page 49
For incompressible (rubber-like) materials, the commonly used models are the neo-Hookean model,
the Mooney-Rivlin model, and the Ogden model, etc [57, 81, 85, 99, 100]. The strain-energy function of
incompressible neo-Hookean materials, for example, is deﬁned by
We
(F) =
µ
2
(tr(FT
F) − 3) or We
(λ1, λ2, λ3) =
µ
2
(λ2
1 + λ2
2 + λ2
3 − 3), (184)
where µ is the shear modulus for small deformation. The constraint of incompressibility reads
det F = 1 or λ1λ2λ3 = 1. (185)
The ﬁrst Piola-Kirchhoﬀ stress Eq. (174) or Eq. (182) is
Te
= µF − LaF−T
or Te
=
3
i=1
µλi −
1
λi
La ˆni ⊗ ˆNi (186)
while the Cauchy stress Eq. (175) or Eq. (183) is
σe
= µFFT
− LaI or σe
=
3
i=1
µλ2
i − La ˆni ⊗ ˆni. (187)
The strain-energy function of the Gent model [101], for example, is
We
(λ1, λ2, λ3) = −
µJm
2
ln 1 −
I1 − 3
Jm
, (188)
where I1 = λ2
1 + λ2
2 + λ2
3 and Jm is a material constant.
4.5.3 Maxwell stress tensor for an ideal dielectric elastomer
In this section, we consider linear dielectrics that admit the relation (149) between the electric ﬁeld E
and the electric displacement D in the current conﬁguration.
Consider a charge density ρf in an electric ﬁeld E. The electric body force is equal to the Lorentz
force ρf E in the electrostatics. By the deﬁnition, div σM
= ρf E, i.e., the divergence of the Maxwell stress
σM
equals to the electric body force ρf E, we have the Maxwell stress in the deformed linear dielectrics
σM
= εE ⊗ E −
ε|E|2
2
I, σM
ij = εEiEj −
ε(E2
1 + E2
2 + E2
3 )
2
δij. (189)
Note that the divergence, div σM
, of the Maxwell stress is equal to the electric body force ρf E in linear
dielectrics.9
9By the identities div (E ⊗ E) = (div E)E + (grad E)E and div (|E|2I) = grad |E|2 = 2(grad E)T E, together with
Eqs. (146)2 and (149), we can obtain the divergence div σM = ρf E + ε[(grad E)E − (grad E)TE]. Since the electric ﬁeld
E in Eq. (147) is irrotational, we have grad E = (grad E)T . Then, (grad E)E = (grad E)T E and the divergence ﬁnally
reduces to div σM = ρf E.
49
Yang & Sharma: A tutorial on the stability and bifurcation Page 50
It follows from Eqs. (165) and (189) that the Piola-Maxwell stress TM
in Eq. (172) is
TM
= εJE ⊗ F−1
E −
εJ|E|2
2
F−T
. (190)
By Eqs. (149) and (154), σM
in Eq. (189) and TM
in Eq. (190) can be, respectively, reformulated as [42]
σM
=
1
ε
D ⊗ D −
|D|2
2ε
I, σM
ij =
1
ε
DiDj −
(D2
1 + D2
2 + D2
3)
2ε
δij, (191)
and
TM
=
1
εJ
F ˜D ⊗ ˜D −
1
2εJ
|F ˜D|2
F−T
. (192)
By the conversion formulas in Table 2, the Piola-Maxwell stress can also be written as either TM
=
1
ε D ⊗ ˜D − J
2ε |D|2
F−T
or TM
= E ⊗ ˜D − εJ
2 |E|2
F−T
.
4.6 Constitutive law for soft dielectrics
The choice of a constitutive model depends on the material being modeled and, of course, the independent
state variables chosen. As an example, in several works, Suo et al. [95] assumed that the free
energy density is a function of the deformation gradient F and the nominal electric displacement ˜D, i.e.,
W(F, ˜D). Thus the nominal stress and the nominal electric ﬁeld are
T =
∂W(F, ˜D)
∂F
and ˜E =
∂W(F, ˜D)
∂ ˜D
. (193)
For the model of ideal soft dielectrics [42,102], the free energy density can be written as
W(F, ˜D) = We
(F) +
1
2εJ
|F ˜D|2
− La(J − 1), (194)
where We
(F) is the strain-energy function due to the mechanical deformation, and La serves as the
Lagrange multiplier related to the constraint of incompressibility, i.e., J = 1. Typically, La is equal to
zero in the case of compressible dielectrics. Substitution of Eq. (194) into Eq. (193) gives
T =
∂W(F, ˜D)
∂F
=
∂We
(F)
∂F
+ TM
− LaF−T
, (195a)
˜E =
∂W(F, ˜D)
∂ ˜D
=
1
εJ
FT
(F ˜D) =
1
εJ
C ˜D, (195b)
where the Maxwell stress TM
is given by Eq. (192). In contrast, the constitutive law can also be
formulated, for example, by using the deformation gradient F and the nominal polarization ˜P, the
deformation gradient F and the nominal electric ﬁeld ˜E.
50
Yang & Sharma: A tutorial on the stability and bifurcation Page 51
4.7 Summary of the boundary-value problem
The 3D BVP equations of electrostatics of deformable dielectrics in the reference conﬁguration can be
summarized as:
Div T + b0R
= 0
equilibrium equation
, Curl ˜E = 0, Div ˜D = ˜ρf
Maxwell’s equations
, ˜D = εJC−1 ˜E
relation (linear dielectrics)
,
u = u0 on ∂Ωu
R, TN = tR on ∂Ωt
R
mechanical BCs
, ξ = ξ0 on ∂Ωξ
R, ˜D · N = Q on ∂Ωd
R
electric BCs
,
T = Te
+ TM
, F = χ(X) = u + I
kinematics
,
Te
=



∂We
(F)
∂F
, compressible solids,
∂We
(F)
∂F
− LaF−T
, J = 1, incompressible solids,
TM
=
1
εJ
F ˜D ⊗ ˜D −
1
2εJ
|F ˜D|2
F−T
Maxwell stress
.



(196)
In the above representation, we have chosen the constitutive law representative of a linear dielectric,
and the nominal electric ﬁeld is our state variable (in addition to deformation).
4.8 Incremental formulation and bifurcation analysis
In analogy with Sec. 3.6, we now linearize the 3D BVP (196) equations. The inﬁnitesimal increment of
the deformation χ is denoted by u (X), i.e., χ(X) → χ(X)+u (X), with suﬃciently small u (X) and
u (X) . In contrast to Eq. (62), the inﬁnitesimal increment of the deformation gradient is
F =
d
dτ
χ(X) + τu (X)
τ=0
= u . (197)
Similarly, it is easy to show that (FT
) = (F )T
. By the variation F → F + F , F 1, we can
51
Yang & Sharma: A tutorial on the stability and bifurcation Page 52
also obtain the following inﬁnitesimal increments10
(F−1
) = −F−1
F F−1
, (F−T
) = −F−T
(F )T
F−T
= [(F−1
) ]T
,
C = FT
F + (FT
) F, (C−1
) = −C−1
C C−1
,
J = JF−T
: F , (J−1
) = −J−2
J .



(198)
The inﬁnitesimal increment of the electric potential ξ is denoted by ξ (X), i.e., ξ(X) → ξ(X)+ξ (X),
with suﬃciently small ξ (X) and ξ (X) . The inﬁnitesimal increment of the nominal electric ﬁeld
˜E is expressed as
˜E =
d
dτ
− ξ(X) + τξ (X)
τ=0
= − ξ . (199)
Consider the second row with variables (˜E, F) in Table 2. The inﬁnitesimal increment ˜D of the
nominal electric displacement ˜D = εJC−1 ˜E is
˜D = εJ C−1 ˜E + εJ(C−1
) ˜E + εJC−1 ˜E , (200)
and the inﬁnitesimal increment E of the true electric ﬁeld E = F−T ˜E is
E = (F−T
) ˜E + F−T ˜E , (201)
where the inﬁnitesimal increments J , (C−1
) , (F−T
) , and ˜E are given in Eqs. (198) and (199).
The detailed derivations are as follows: ˜D = d
dτ ε(J + τJ )(C−1
+ τ(C−1
) )(˜E + τ ˜E )
τ=0
and
E = d
dτ (F−T
+ τ(F−T
) )(˜E + τ ˜E )
τ=0
.
The inﬁnitesimal increment of the Maxwell stress, i.e., TM
= E ⊗ ˜D − εJ
2 |E|2
F−T
, then is
(TM
) = E ⊗ ˜D + E ⊗ ˜D −
ε
2
J |E|2
F−T
+ 2J(E · E )F−T
+ J|E|2
(F−T
) . (202)
10By the deﬁnition (13), the detailed derivations of the inﬁnitesimal increments are listed as follows:
(F−1) = d
dτ
(F + τF )−1
τ=0
= limτ→0
(F+τF )−1
−F−1
τ
= limτ→0
(F+τF )−1
{F−(F+τF )}F−1
τ
=
−F−1F F−1, (F−T ) = limτ→0
(F+τF )−T
−F−T
τ
= limτ→0
(F+τF )−T
FT
−(F+τF )T
F−T
τ
=
−F−T (F )T F−T , C = (FT F) = d
dτ
(F + τF )T
(F + τF )
τ=0
= limτ→0
(F+τF )T
(F+τF )−FT
F
τ
=
FT F + (F )T F = FT F + (FT ) F, (C−1) = (F−1F−T ) = limτ→0
(F+τF )−1
(F+τF )−T
−F−1
F−T
τ
=
limτ→0
(F+τF )−1
(F+τF )−T
FT
F−(F+τF )T
(F+τF ) F−1
F−T
τ
= −C−1C C−1, J = (det F) =
limτ→0
det(F+τF )−det F
τ
= JF−T : F . Here we use the relation det (F + τF ) = (det F) det I + τF−1F = J[1 +
τ(F−T : F )+o(τ)], which comes from the Cayley-Hamilton equation, cf., e.g., page 35 in [80]. From the forms of the above
increments, we have C ∼ F 1, (F−1) ∼ F 1, (C−1) ∼ C ∼ F 1, and J ∼ F 1.
By the variation J → J + J , J 1, we have (J−1) = limτ→0
(J+τJ )−1
−J−1
τ
= limτ→0
J−(J+τJ )
τJ(J+τJ )
= −J−2J ,
and, hence, that (J−1) ∼ J ∼ F 1.
52
Yang & Sharma: A tutorial on the stability and bifurcation Page 53
Now the incremental version of the boundary-value problem (196) can be represented by
Div T = 0, Curl ˜E = 0, Div ˜D = 0,
u = 0 on ∂Ωu
R, T N = 0 on ∂Ωt
R
incremental mechanical BCs
, ξ = 0 on ∂Ωξ
R, ˜D · N = 0 on ∂Ωd
R
incremental electric BCs
,
T = (Te
) + (TM
) ,
(Te
) =



∂2
We
(F)
∂F2
F , compressible solids,
∂2
We
(F)
∂F2
F − LaF−T
− La(F−T
) , J = 0 for incompressible solids,



(203)
along with Eqs. (197)−(202).
4.9 Example 1: Pull-in instability of a dielectric elastomer ﬁlm by using the
incremental method
Figure 16: Schematic of a deformed ﬁlm of dielectric elastomers coated with two compliant electrodes.
The dead loads P1 and P2, and the voltage Φ deform the ﬁlm from L1, L2, and L3 to λ1L1, λ2L2, and
λ3L3, as well as induce an electric charge of magnitude Q on either electrode [38].
Consider an undeformed dielectric elastomer ﬁlm with the dimensions L1, L2, and L3. The upper
and bottom surfaces of the ﬁlm are both coated with compliant electrodes. The ﬁlm deforms from
L1 × L2 × L3 to λ1L1 × λ2L2 × λ3L3 by applying the in-plane dead loads P1 and P2 and a voltage Φ
across its thickness (see Fig. 16).
Taking the Cartesian coordinates (X1, X2, X3) with an orthonormal basis (e1, e2, e3), the domain
53
Yang & Sharma: A tutorial on the stability and bifurcation Page 54
occupied by the undeformed ﬁlm is represented by
ΩR = (X1, X2, X3) ∈ R3
: 0 ≤ X1 ≤ L1, 0 ≤ X2 ≤ L2, 0 ≤ X3 ≤ L3 . (204)
The boundary ∂ΩR of ΩR consists of six parts:
S1− = X ∈ R3
: X1 = 0 , S1+ = X ∈ R3
: X1 = L1 ,
S2− = X ∈ R3
: X2 = 0 , S2+ = X ∈ R3
: X2 = L2 ,
S3− = X ∈ R3
: X3 = 0 , S3+ = X ∈ R3
: X3 = L3 .



(205)
The mechanical and electric boundary conditions are
T(−e1) = tR = −s1e1, ˜D · (−e1) = 0 on S1−, Te1 = tR = s1e1, ˜D · e1 = 0 on S1+,
T(−e2) = tR = −s2e2, ˜D · (−e2) = 0 on S2−, Te2 = tR = s2e2, ˜D · e2 = 0 on S2+,
T(−e3) = tR = 0, ξ = Φ on S3−, Te3 = tR = 0, ξ = 0 on S3+,



(206)
where s1 = P1/(L2L3) and s2 = P2/(L1L3).
4.9.1 Large deformation
We assume that the dielectric ﬁlm admits the following homogeneous deformation
x1 = λ1X1, x2 = λ2X2, x3 = λ3X3, (207)
where λ1, λ2, and λ3 are constant stretches. Then the matrix representation of the deformation gradient
is a 3 × 3 diagonal matrix
F :=





λ1 0 0
0 λ2 0
0 0 λ3





= diag (λ1, λ2, λ3). (208)
Consequently,
F−1
= F−T
:= diag (λ−1
1 , λ−1
2 , λ−1
3 ), C−1
:= diag (λ−2
1 , λ−2
2 , λ−2
3 ), J = det F = λ1λ2λ3 = 1. (209)
Note that F, F−1
, F−T
and C−1
are independent of the coordinates. For linear dielectric elastomers
and in absence of free changes, the Maxwell equations in the BVP (196) yield the three-dimensional
Laplace equation
∂2
ξ
∂X2
1
+
∂2
ξ
∂X2
2
+
∂2
ξ
∂X2
3
= 0. (210)
By the electric boundary conditions in Eq. (206), the solution of ξ is given by
ξ(X) = Φ
L3 − X3
L3
, (211)
54
Yang & Sharma: A tutorial on the stability and bifurcation Page 55
and the nominal electric ﬁeld and displacement are
˜E = − ξ =
Φ
L3
e3 = ˜Ee3 and ˜D = εJC−1 ˜E = ελ−2
3
˜Ee3 = ˜De3 (212)
that are independent of the coordinates. By Eqs. (208) and (211), the Maxwell stress is obtained as
TM
:=
ε ˜E2
2
diag (−λ−1
1 λ−2
3 , −λ−1
2 λ−2
3 , λ−3
3 ). (213)
For incompressible neo-Hookean materials, the elastic stress is given by Eq. (186). By Eqs. (208),
(209) and (213), the total stress tensor is
T := diag (µλ1 − Laλ−1
1 −
1
2
ελ−1
1 λ−2
3
˜E2
, µλ2 − Laλ−1
2 −
1
2
ελ−1
2 λ−2
3
˜E2
, µλ3 − Laλ−1
3 +
1
2
ελ−3
3
˜E2
). (214)
By the mechanical boundary conditions in Eq. (206), we ﬁnally have
µλ1 − Laλ−1
1 −
1
2
ελ−1
1 λ−2
3
˜E2
= s1, (215a)
µλ2 − Laλ−1
2 −
1
2
ελ−1
2 λ−2
3
˜E2
= s2, (215b)
µλ3 − Laλ−1
3 +
1
2
ελ−3
3
˜E2
= 0, (215c)
which, together with the constraint of incompressibility, i.e., λ1λ2λ3 = 1, can provide the solutions
of (λ1, λ2, λ3, La). We can also reduce the number of the algebraic equations from four to two. By
substituting λ3 = λ−1
1 λ−1
2 and La = µλ2
3 + 1
2 ελ−2
3
˜E2
into Eqs. (215a) and (215a), we have
µ(λ1 − λ−3
1 λ−2
2 ) − ελ1λ2
2
˜E2
= s1, (216a)
µ(λ2 − λ−3
2 λ−2
1 ) − ελ2λ2
1
˜E2
= s2, (216b)
which, by the relation ˜E = ε−1
λ−2
1 λ−2
2
˜D in Eq. (212)2, are the same as Eqs. (6a) and (6b) in [38].
A trivial solution to the boundary-value problem is summarized here: the nominal electric ﬁeld
is ˜E = Φ
L3
e3 = ˜Ee3 in Eq. (212), stretches λ1 and λ2 are the solutions of Eqs. (216a) and (216b),
λ3 = λ−1
1 λ−1
2 and the Lagrange multiplier is La = µλ2
3 + 1
2 ελ−2
3
˜E2
.
4.9.2 Incremental boundary-value problem
In contrast to Eq. (206), the incremental mechanical and electric boundary conditions are
T (−e1) = 0, ˜D · (−e1) = 0 on S1−, T e1 = 0, ˜D · e1 = 0 on S1+,
T (−e2) = 0, ˜D · (−e2) = 0 on S2−, T e2 = 0, ˜D · e2 = 0 on S2+,
T (−e3) = 0, ξ = 0 on S3−, T e3 = 0, ξ = 0 on S3+.



(217)
55
Yang & Sharma: A tutorial on the stability and bifurcation Page 56
For the trivial solution, the inﬁnitesimal increments of the displacement are assumed to be
u1 = λ1X1, u2 = λ2X2, u3 = λ3X3, (218)
where λ1, λ2, and λ3 are suﬃciently small and independent of the coordinates. By Eqs. (197) and (218),
the inﬁnitesimal increment of the deformation gradient is
F := diag (λ1, λ2, λ3). (219)
With Eqs. (208) and (219), other increments in Eq. (198) for the trivial solution can be obtained as
(F−1
) = (F−T
) := −diag (λ−2
1 λ1, λ−2
2 λ2, λ−2
3 λ3), J = λ2λ3λ1 + λ3λ1λ2 + λ1λ2λ3,
C := 2 diag (λ1λ1, λ2λ2, λ3λ3), (C−1
) := −2 diag (λ−3
1 λ1, λ−3
2 λ2, λ−3
3 λ3).



(220)
For incompressible solids, the inﬁnitesimal increment of the Jacobian is zero, i.e.,
J = λ2λ3λ1 + λ3λ1λ2 + λ1λ2λ3 = 0. (221)
Consider incompressible neo-Hookean solids (184). By Eqs. (208) and (220), the incremental elastic
stress tensor in the incremental BVP (203) is obtained as
(Te
) := diag (µλ1 − λ−1
1 La + Laλ−2
1 λ1, µλ2 − λ−1
2 La + Laλ−2
2 λ2, µλ3 − λ−1
3 La + Laλ−2
3 λ3), (222)
where La = µλ2
3 + 1
2 ελ−2
3
˜E2
and La is the inﬁnitesimal increment of the Lagrange multiplier.
The inﬁnitesimal increment of the electric potential is ξ , then the incremental nominal electric ﬁeld
(199) is ˜E = − ξ . Also, the incremental nominal electric displacement (200) is ˜D = ε(C−1
) ˜E +
εC−1 ˜E by using J = 1 and J = 0. Since C−1
in Eq. (209), (C−1
) in Eq. (220), ˜E in Eq. (212)
are independent of the coordinates, the incremental Maxwell equation Div ˜D = 0 in Eq. (203) implies
Div ˜E = 0, namely
∂2
ξ
∂X2
1
+
∂2
ξ
∂X2
2
+
∂2
ξ
∂X2
3
= 0, (223)
which, together with the incremental electric boundary conditions in Eq. (217), gives the solution
ξ (X) = 0. (224)
Then the increments from the trivial solution become
˜E = 0 and ˜D = −2ελ−3
3
˜Eλ3e3. (225)
By Eqs. (209)1 and (212)1, the true electric ﬁeld is E = F−T ˜E = λ−1
3
˜Ee3. It follows from Eqs. (220)1,
(212)1 and (225)1 that the increment E in Eq. (201) is obtained as E = −λ−2
3
˜Eλ3e3. Thus, the
56
Yang & Sharma: A tutorial on the stability and bifurcation Page 57
incremental Maxwell stress (202) is
(TM
) :=
ε ˜E2
2
diag (λ−2
1 λ−2
3 λ1 + 2λ−1
1 λ−3
3 λ3, λ−2
2 λ−2
3 λ2 + 2λ−1
2 λ−3
3 λ3, −3λ−4
3 λ3). (226)
By Eqs. (222) and (226), we have the increment T of the total stress in the incremental BVP (203).
Thus, the incremental mechanical boundary conditions in Eq. (217) give
(µ + µλ−2
1 λ2
3 + ελ−2
1 λ−2
3
˜E2
)λ1 + ελ−1
1 λ−3
3
˜E2
λ3 − λ−1
1 La = 0, (227a)
(µ + µλ−2
2 λ2
3 + ελ−2
2 λ−2
3
˜E2
)λ2 + ελ−1
2 λ−3
3
˜E2
λ3 − λ−1
2 La = 0, (227b)
(2µ − ελ−4
3
˜E2
)λ3 − λ−1
3 La = 0, (227c)
which, together with Eq. (221), can give the solutions of (λ1, λ2, λ3, La).
4.9.3 Solution of the incremental boundary-value problem
Here we have 4 equations in 4 unknowns, i.e., Eqs. (221), (227a)−(227c) with variables (λ1, λ2, λ3, La). It
follows from Eq. (227c) that La = (2µλ3 − ελ−3
3
˜E2
)λ3. By substituting La into Eqs. (227a) and (227b),
we can reduce the set of simultaneous equations to 3 equations in 3 unknowns. The set of equations may
be expressed as








µ + µλ−2
1 λ2
3 + ελ−2
1 λ−2
3
˜E2
0 2ελ−1
1 λ−3
3
˜E2
− 2µλ−1
1 λ3
0 µ + µλ−2
2 λ2
3 + ελ−2
2 λ−2
3
˜E2
2ελ−1
2 λ−3
3
˜E2
− 2µλ−1
2 λ3
λ2λ3 λ3λ1 λ1λ2
















λ1
λ2
λ3








= 0.
(228)
If the determinant of the 3 × 3 matrix is non-zero, the simultaneous linear equations only have the
zero solution, i.e., λ1 = λ2 = λ3 = 0. If the determinant is zero, there exist an inﬁnite number of
solutions. Thus, the necessary condition for the existence of non-trivial incremental solutions is a zero
determinant.
In the following, we will show that the condition for a zero determinant of the 3 × 3 matrix in
Eq. (228) is the same as that of the Hessian matrix in [38]. Consider the 3 × 3 matrix in Eq. (228)
as Matrix(I). Taking λ3 = λ−1
1 λ−1
2 and replacing row (1) by row (1) + 2µλ−3
1 λ−2
2 × row (3), row (2)
by row (2) + 2µλ−2
1 λ−3
2 × row (3) in Matrix(I), we obtain Matrix(II). Subsequently, replacing column
(1) by column (1) + λ−2
1 λ−1
2 × column (3), column (2) by column (2) + λ−1
1 λ−2
2 × column (3) in Matrix(II),
we get Matrix(III). We multiply column (3) in Matrix(III) by −ε−1 ˜E−1
λ−3
1 λ−3
2 and obtain
Matrix(IV). Multiplying row (3) in Matrix(IV) by − ˜E, we ﬁnally get Matrix(V), which, by the relation
˜E = ε−1
λ−2
1 λ−2
2
˜D, is exact the Hessian matrix in [38]. By the eﬀects of row/column operations of determinants
[98], we have the following relation: det (Matrix(I)) = det (Matrix(II)) = det (Matrix(III)) =
det (Matrix(IV))/(−ε−1 ˜E−1
λ−3
1 λ−3
2 ) and det (Matrix(IV)) = det (Matrix(V))/(− ˜E). Thus, the condi-
57
Yang & Sharma: A tutorial on the stability and bifurcation Page 58
tion for a zero determinant of the 3 × 3 matrix in Eq. (228) is equivalent to that of the Hessian matrix
in [38].
4.10 Example 2: Wrinkle surface instability of a dielectric elastomer by using
the incremental method
In addition to pull-in instability for a homogeneous deformation, surface instabilities, which lead to inhomogeneous
deformation (wrinkling), are also common. Given the existence of an extensive body of work
on the electromechanical instabilities related to inhomogeneous deformations, we simply point to a small
sample of the works that the reader may consult (and the references therein): [19,22,29,103–118,155].
The schematic of the onset of surface wrinkling induced by an applied electric ﬁeld is shown in Fig. 17.
Due to inhomogeneous deformations, the electro-wrinkling problem is more involved. We will use linear
bifurcation analysis to ﬁnd the necessary condition for the existence of a non-trivial solution (surface
wrinkling). We employ a similar procedure as what can be found in the appendix of the work by Wang
and Zhao [19].
(a) (b)
Figure 17: Schematic of surface wrinkling of a dielectric block bonded to a rigid substrate. The upper
surface is bonded with a compliant electrode while the rigid substrate is coated with a metal electrode. A
voltage diﬀerence is applied across the two electrodes. This example is inspired by [19]. (a) Undeformed
state with a ﬂat surface. (b) Wrinkled state by increasing the voltage at the threshold.
4.10.1 Large deformation
For the electromechanical system in Fig. 17, the displacement of the bottom surface is ﬁxed. We consider
a 2D formulation within the plane strain assumption. Thus the kinematical boundary conditions at the
bottom surface are:
u(X) = 0 or u1(X1, X2) = u2(X1, X2) = 0 at X2 = 0. (229)
58
Yang & Sharma: A tutorial on the stability and bifurcation Page 59
The deformation gradient is
F = (X + u) = I + u :=


1 + u1,1 u1,2
u2,1 1 + u2,2

 . (230)
The constraint of incompressibility reads
J = det F =
1 + u1,1 u1,2
u2,1 1 + u2,2
= 1. (231)
Alternatively, calculation of the determinant (231) gives
u1,1 + u2,2 + u1,1u2,2 − u1,2u2,1 = 0. (232)
The traction-free boundary conditions at the upper surface are11
Te2 = 0 or T12(X1, X2) = T22(X1, X2) = 0 at X2 = H. (233)
The true electric ﬁeld induced by an applied voltage Φ can be approximately written as
E = Eδ e2 = (0, Eδ)T
, (234)
where
Eδ =
Φ
H + δH
=
Φ
H + u2(X1, H)
. (235)
The total Cauchy stress σ in Eq. (173) is rewritten here as
σ = σe
+ σM
, (236)
where the elastic stress and the Maxwell stress for neo-Hookean ideal dielectrics are
σe
= µFFT
− LaI = µ(I + u + uT
+ u uT
) − LaI, (237)
and
σM
= εE ⊗ E −
ε|E|2
2
I :=
ε
2
E2
δ


−1 0
0 1

 . (238)
In contrast to Eq. (236), the (ﬁrst) Piola-Kirchhoﬀ stress, with the relation (170), is given by
T = JσF−T
= (σe
+ σM
)F−T
, (239)
11Te2 = (T11e1 ⊗ e1 + T12e1 ⊗ e2 + T21e2 ⊗ e1 + T22e2 ⊗ e2)e2 = T12e1 + T22e2 = 0.
59
Yang & Sharma: A tutorial on the stability and bifurcation Page 60
where J = 1 for incompressible materials. Without the body force, the equilibrium equation is
Div T = 0. (240)
Now we have constructed a boundary-value problem:
Div T = 0, det F = 1, E = Eδe2,
u1(X1, X2) = u2(X1, X2) = 0 at X2 = 0,
T12(X1, X2) = T22(X1, X2) = 0 at X2 = H,



(241)
with
F = I + u, T = (σe
+ σM
)F−T
, Eδ =
Φ
H + u2(X1, H)
,
σe
= µ(I + u + uT
+ u uT
) − LaI, σM
:=
ε
2
E2
δ


−1 0
0 1

 .



(242)
A trivial solution to the boundary-value problem is
u0 = 0, La0 = µ +
ε
2
Φ
H
2
, (243)
implying
F0 = I, σe
0 :=
ε
2
Φ
H
2


−1 0
0 −1

 , σM
0 :=
ε
2
Φ
H
2


−1 0
0 1

 . (244)
4.10.2 Small deformation simpliﬁcation
Prior to subsequent stability analysis, it is instructive to ﬁrst simplify the problem for small deformation.
The constraint of incompressibility (232) reduces to
u1,1 + u2,2 = 0 or Div u = 0, (245)
the magnitude of the true electric ﬁeld (235) becomes
Eδ =
Φ
H
1 −
u2(x1, H)
H
, (246)
the elastic stress (237) and the Maxwell stress (238) are reduced to
σe
= µ(I + u + uT
) − LaI, (247)
σM
:=
εE2
0
2
1 − 2
u2(x1, H)
H


−1 0
0 1

 . (248)
60
Yang & Sharma: A tutorial on the stability and bifurcation Page 61
Finally, the boundary-value problem (241) and (242) reduces to
Div σ = 0, Div u = 0, Eδ = E0 1 −
u2(x1, H)
H
,
u1(x1, x2) = u2(x1, x2) = 0 at x2 = 0,
σ12(x1, x2) = σ22(x1, x2) = 0 at x2 = H,



(249)
with
σ = σe
+ σM
, σe
= µ(I + u + uT
) − LaI,
σM
:=
εE2
0
2
1 − 2
u2(x1, H)
H


−1 0
0 1

 .



(250)
4.10.3 Incremental boundary-value problem
For the trivial solution (243), the incremental forms of the displacement and the Lagrange multiplier are
u = u0 + u = u , La = La0 + La, (251)
where u and La are increments with suﬃciently small u , u , La , and La . Then the
incremental form of the constraint Div u = 0 in Eq. (249) is
Div u = u1,1 + u2,2 = 0. (252)
The incremental form of the elastic stress σe
in Eq. (250) is
σe
= σe
0 + (σe
) = σe
0 + µ( u + ( u )T
) − LaI. (253)
The magnitude of the electric ﬁeld in Eq. (249) is
Eδ = E0 1 −
u2(x1, H)
H
(254)
and the incremental form of the Maxwell stress is
σM
= σM
0 + (σM
) :=
ε
2
E2
0


−1 0
0 1

 − εE2
0
u2(x1, H)
H


−1 0
0 1

 . (255)
Therefore, the incremental form of the boundary-value problem (249) and (250) is
Div σ = 0, u1,1 + u2,2 = 0, Eδ = −E0
u2(x1, H)
H
,
u1(x1, x2) = u2(x1, x2) = 0 at x2 = 0,
σ12(x1, x2) = σ22(x1, x2) = 0 at x2 = H,



(256)
61
Yang & Sharma: A tutorial on the stability and bifurcation Page 62
with
σ = (σe
) + (σM
) , (σe
) = µ( u + ( uT
) ) − LaI,
(σM
) := −εE2
0
u2(x1, H)
H


−1 0
0 1

 .



(257)
4.10.4 Solution of the incremental boundary-value problem
An approximation used here is that Div (σM
) = 0 [19, 108]. With the constraint, u1,1 + u2,2 = 0, we
have Div ( uT
) = 0.12
Thus, the equilibrium equation in Eq. (256) reduces to
Div (µ u − LaI) = µ 2
u − La = 0. (258)
With the constraint u1,1 + u2,2 = 0, the general solution of u (x1, x2) can be written as
u1(x1, x2) =
∂φ
∂x2
and u2(x1, x2) = −
∂φ
∂x1
, (259)
where φ(x1, x2) is an arbitrary function. Substitution of Eq. (259) into Eq. (257) gives
σ :=







2µ
∂2
φ
∂x2∂x1
− La µ
∂2
φ
∂x2
2
−
∂2
φ
∂x2
1
µ
∂2
φ
∂x2
2
−
∂2
φ
∂x2
1
−2µ
∂2
φ
∂x1∂x2
− La







+
εE2
0
H
∂φ(x1, H)
∂x1


−1 0
0 1

 . (260)
By using the general solution (259), the equilibrium equation (258) becomes
µ
∂3
φ
∂x2∂x2
1
+
∂3
φ
∂x3
2
−
∂La
∂x1
= 0, µ
∂3
φ
∂x3
1
+
∂3
φ
∂x1∂x2
2
+
∂La
∂x2
= 0, (261)
the kinematical boundary conditions in Eq. (256) read
∂φ
∂x2
=
∂φ
∂x1
= 0 at x2 = 0, (262)
and the traction-free boundary conditions in Eq. (256) are
∂2
φ
∂x2
2
−
∂2
φ
∂x2
1
= 0, −2µ
∂2
φ
∂x1∂x2
− La +
εE2
0
H
∂φ(x1, H)
∂x1
= 0 at x2 = H. (263)
The boundary-value problem (261)−(263) is quite classical and the reader can refer to the following
works for more details: [19,35,108,119]. Consider the forms of the solution
φ(x1, x2) = ¯φ(x2) cos(kx1) and La(x1, x2) = ¯L(x2) sin(kx1). (264)
12Since u1,1 + u2,2 = ui,i = uj,j = 0, i, j = 1, 2, then Div ( uT ) = Div (uj,iei ⊗ ej) = uj,ijei = ∂
∂xi
uj,jei = 0.
62
Yang & Sharma: A tutorial on the stability and bifurcation Page 63
Then the boundary-value problem (261)−(263) is converted to
µ − k2 ∂ ¯φ
∂x2
+
∂3 ¯φ
∂x3
2
− k ¯L = 0, µ k3 ¯φ − k
∂2 ¯φ
∂x2
2
+
∂ ¯L
∂x2
= 0, x2 ∈ (0, H), (265a)
∂ ¯φ
∂x2
= 0, ¯φ = 0 at x2 = 0, (265b)
∂2 ¯φ
∂x2
2
+ k2 ¯φ = 0, 2µk
∂ ¯φ
∂x2
− ¯L −
εE2
0
H
k ¯φ = 0 at x2 = H. (265c)
By eliminating ¯L in Eq. (265a), we have a fourth-order diﬀerential equation
∂4 ¯φ
∂x4
2
− 2k2 ∂2 ¯φ
∂x2
2
+ k4 ¯φ = 0 (266)
with the general solution
¯φ(x2) = (C1 + C2x2)ekx2
+ (C3 + C4x2)e−kx2
, (267)
where Ci, i = 1, 2, 3, 4, are constant. Substituting Eq. (267) into Eq. (265a)1, we have
¯L(x2) = 2µk(C2ekx2
+ C4e−kx2
). (268)
Consider the boundary conditions (265b) and (265c). We encounter a set of simultaneous linear
equations, that is, we have four equations in four unknowns C1, C2, C3, C4 of the form








k 1 −k 1
1 0 1 0
2k2
ekH
2(1 + kH)kekH
2k2
e−kH
2(−1 + kH)ke−kH
(
εE2
0
H
− 2kµ)kekH
(εE2
0 − 2kHµ)kekH
(
εE2
0
H
+ 2kµ)ke−kH
(εE2
0 + 2kHµ)ke−kH
















C1
C2
C3
C4








= 0.
(269)
The above equations (269) only have non-trivial solutions Ci, i = 1, 2, 3, 4, if the determinant of the
coeﬃcient matrix vanishes, which gives
εE2
0
µ
=
2kH e4kH
+ (2 + 4k2
H2
)e2kH
+ 1
e4kH − 4kHe2kH − 1
. (270)
Equation (270) is exactly equation (1) in the work by Wang and Zhao [19] in the absence of surface
energy. In Fig. 18, we plot how the electric ﬁeld E0 ε/µ varies with respect to the ratio kH. The lowest
electric ﬁeld (the threshold) for the onset of electromechanical surface wrinkling is E0 ε/µ = 2.49 at
kH = 2.12. If the applied electric ﬁeld is less than the threshold, there is no electromechanical wrinkling
and the elastomer has a ﬂat surface.
The block dimensions would indeed aﬀect the critical electric ﬁeld and the wavelength 2π/k (or the
wavenumber k) of surface wrinkling. In the analysis, we regard a continuous wavenumber k to carry out
63
Yang & Sharma: A tutorial on the stability and bifurcation Page 64
Figure 18: The electric ﬁeld for the onset of surface wrinkling of neo-Hookean ideal dielectrics. The
(blue) curve comes from the numerical plot of Eq. (270). The lowest point on the curve is denoted by a
cross mark at E0 ε/µ = 2.49 and kH = 2.12. The cross mark corresponds to the lowest electric ﬁeld
and can be regarded as the threshold of the surface wrinkling instability. Below the lowest point of the
curve, the elastomer has a ﬂat surface, which means that the system shown in Fig. 17 is stable and there
is no electromechanical wrinkling.
our plot. This treatment is accurate enough if the length of the elastomer block is suﬃciently large, which
admits relatively dense discrete values of k. In addition to the length, the height H of the elastomer
block aﬀects the wrinkling phenomenon in terms of the product kH. For the onset of surface wrinkling
in Eq. (270), a larger electric ﬁeld corresponds to a smaller product kH, giving a larger wavelength 2π/k
of wrinkles for an elastomer block with a given geometry. A detailed discussion of the eﬀects of block
dimensions and side boundary conditions on the surface wrinkling is beyond the scope of this tutorial.
5 Electromechanical instability and bifurcation of soft dielectrics:
pull-in instability
In this section, we mainly revisit pull-in instability of soft dielectrics by using the energy method. For
other types of instability, the reader may refer to the following two excellent reviews [42,59]. There is an
extensive body of work on pull-in instability of a dielectric elastomer ﬁlm for homogeneous deformation
[38,40,120–124]. Here we focus on two works [38,40] and discuss pull-in instability by both the principle
of minimum energy (i.e., the energy method) and using linear bifurcation analysis.
64
Yang & Sharma: A tutorial on the stability and bifurcation Page 65
5.1 Pull-in instability of a dielectric elastomer ﬁlm
We review the example of electromechanical instability of a ﬁlm of a dielectric elastomer by Zhao and
Suo [38]. The stability is examined by using the principle of minimum energy. Since, in this example,
the deformation is homogeneous, this is the same as the so-called Hessian approach.
5.1.1 Stability analysis by using the Hessian approach
Consider the electrostatic system shown in Fig. 16. Subject to in-plane dead loads P1 and P2 and
in-thickness voltage Φ, the ﬁlm deforms from L1 × L2 × L3 to λ1L1 × λ2L2 × λ3L3. The constraint
of incompressibility gives the relation λ3 = 1/(λ1λ2) between the stretches. The free energy G of the
system is given by [38,95]
G
L1L2L3
= W(λ1, λ2, ˜D) −
P1
L2L3
λ1 −
P2
L1L3
λ2 −
Φ
L3
˜D, (271)
where W(λ1, λ2, ˜D) is the free energy function. With the small variations of the generalized coordinates,
δλ1, δλ2, δ ˜D, the variation of the free energy is
∆G
L1L2L3
=
∂W
∂λ1
−
P1
L2L3
δλ1 +
∂W
∂λ2
−
P2
L1L3
δλ2 +
∂W
∂ ˜D
−
Φ
L3
δ ˜D
+
1
2
∂2
W
∂λ2
1
(δλ1)2
+
1
2
∂2
W
∂λ2
2
(δλ2)2
+
1
2
∂2
W
∂ ˜D2
(δ ˜D)2
+
∂2
W
∂λ1∂λ2
δλ1δλ2 +
∂2
W
∂λ1∂ ˜D
δλ1δ ˜D +
∂2
W
∂λ2∂ ˜D
δλ2δ ˜D, (272)
where only the ﬁrst and second variations are retained and all the high-order terms are omitted.
In equilibrium, the ﬁrst variation is zero, we have
∂W
∂λ1
−
P1
L2L3
= 0,
∂W
∂λ2
−
P2
L1L3
= 0,
∂W
∂ ˜D
−
Φ
L3
= 0. (273)
In Eq. (272), the second variation is related through the real-symmetric Hessian matrix
H =










∂2
W
∂λ2
1
∂2
W
∂λ1∂λ2
∂2
W
∂λ1∂ ˜D
∂2
W
∂λ2∂λ1
∂2
W
∂λ2
2
∂2
W
∂λ2∂ ˜D
∂2
W
∂ ˜D∂λ1
∂2
W
∂ ˜D∂λ2
∂2
W
∂ ˜D2










, (274)
which must be positive deﬁnite if the equilibrium state is stable. A real-symmetric positive deﬁnite
matrix requires that all the principal minors of the matrix are positive [98]. The 1 × 1 principal minors
are the diagonal entries
∂2
W
∂λ2
1
,
∂2
W
∂λ2
2
,
∂2
W
∂ ˜D2
. (275)
The 2 × 2 principal minors are
65
Yang & Sharma: A tutorial on the stability and bifurcation Page 66
∂2
W
∂λ2
1
∂2
W
∂λ1∂λ2
∂2
W
∂λ2∂λ1
∂2
W
∂λ2
2
,
∂2
W
∂λ2
1
∂2
W
∂λ1∂ ˜D
∂2
W
∂ ˜D∂λ1
∂2
W
∂ ˜D2
,
∂2
W
∂λ2
2
∂2
W
∂λ2∂ ˜D
∂2
W
∂ ˜D∂λ2
∂2
W
∂ ˜D2
, (276)
and the only 3 × 3 principal minor is
det H =
∂2
W
∂λ2
1
∂2
W
∂λ1∂λ2
∂2
W
∂λ1∂ ˜D
∂2
W
∂λ2∂λ1
∂2
W
∂λ2
2
∂2
W
∂λ2∂ ˜D
∂2
W
∂ ˜D∂λ1
∂2
W
∂ ˜D∂λ2
∂2
W
∂ ˜D2
. (277)
For an ideal dielectric elastomer (194), the free energy function is given by
W(λ1, λ2, ˜D) =
µ
2
(λ2
1 + λ2
2 + λ−2
1 λ−2
2 − 3) +
˜D2
2ε
λ−2
1 λ−2
2 . (278)
Substitution of Eq. (278) into the equilibrium equations (273) and the Hessian matrix (274) yields
the necessary equations for the stability analysis as can be found in the works [38,121].
For notational expediency, we summarize the stability analysis in vector form. The free energy in
Eq. (271) and the variation in Eq. (272) can be written as
G(x; L) and G(x + δx; L) = G(x; L) +
∂G
∂x
· δx + δx ·
∂2
G
∂x2
[δx] + o(δx2
), (279)
where the state is x = (λ1, λ2, ˜D)T
, the increments are δx = (δλ1, δλ2, δ ˜D)T
, and the load parameters
are L = (P1, P2, Φ)T
. The function o(δx2
) has the limiting behavior, i.e., the ratio o(δx2
)
δx2 approaches
zero as δx2
→ 0.
The stability of the state x subjected to the load parameter L requires a local energy minimization,
such that G(x; L) ≤ G(x + δx; L) for suﬃciently small δx . Thus we have the equilibrium equations
∂G(x; L)
∂x
· δx = 0 =⇒ gj(x; L) =
∂G(x; L)
∂xj
= 0, (280)
which correspond to Eq. (273), and the stability conditions
δx ·
∂2
G
∂x2
[δx] ≥ 0 =⇒



a Hessian matrix H =
∂2
G
∂x2
is (semi) positive deﬁnite
and the matrix elements are Hij =
∂2
G(x; L)
∂xi∂xj
,
(281)
where i, j, k = 1, 2, 3.
66
Yang & Sharma: A tutorial on the stability and bifurcation Page 67
It follows from Eqs. (280), (278), and (271) that the equilibrium equations are
µ(λ1 − λ−3
1 λ−2
2 ) −
˜D2
ε
λ−3
1 λ−2
2 −
P1
L2L3
= 0, (282a)
µ(λ2 − λ−2
1 λ−3
2 ) −
˜D2
ε
λ−2
1 λ−3
2 −
P2
L1L3
= 0, (282b)
˜D
ε
λ−2
1 λ−2
2 −
Φ
L3
= 0. (282c)
It follows from Eqs. (281), (278), and (271) that the Hessian matrix is
H =









µ(1 + 3λ−4
1 λ−2
2 ) +
3 ˜D2
ε
λ−4
1 λ−2
2 2µλ−3
1 λ−3
2 +
2 ˜D2
ε
λ−3
1 λ−3
2 −
2 ˜D
ε
λ−3
1 λ−2
2
2µλ−3
1 λ−3
2 +
2 ˜D2
ε
λ−3
1 λ−3
2 µ(1 + 3λ−2
1 λ−4
2 ) +
3 ˜D2
ε
λ−2
1 λ−4
2 −
2 ˜D
ε
λ−2
1 λ−3
2
−
2 ˜D
ε
λ−3
1 λ−2
2 −
2 ˜D
ε
λ−2
1 λ−3
2
1
ε
λ−2
1 λ−2
2









. (283)
As stated in Eq. (281), the Hessian matrix (283) is positive deﬁnite when the electrostatic system
shown in Fig. 16 at equilibrium (282) is stable.
5.1.2 Bifurcation analysis
We now discuss the bifurcation problem pertaining to the electromechanical system shown in Fig. 16.
Based on the equilibrium equations (273), we will formulate the incremental equilibrium equations and
then ﬁnd the incremental solutions.13
The construction of the incremental equilibrium equations can be
found in the following books [57,58].
Let us describe the idea of linear bifurcation analysis in general terms ﬁrst and then we discuss the
bifurcation of the system shown in Fig. 16 in detail. Consider the equilibrium equations (280). A small
perturbation of the equilibrium state x = (λ1, λ2, ˜D)T
is given by δx = (δλ1, δλ2, δ ˜D)T
for the load
parameters L = (P1, P2, Φ)T
. The new state is x + δx = (λ1 + δλ1, λ2 + δλ2, ˜D + δ ˜D)T
with suﬃciently
small δx . Then the expansion of the function gj in Eq. (280) for the new state is
gj(x + δx; L) = gj(x; L) +
∂gj
∂x (x;L)
· δx + o(δx). (284)
By the implicit function theorem in Sec. 2.4.1, a necessary condition for (x; L) to be a bifurcation
point is that the Fr´echet derivative
∂gj
∂x (x;L)
(285)
not be invertible. Otherwise, there exists a function h, such that gj(h(L); L) = 0.
13The method of incremental solutions is based on the adjacent equilibrium stability criterion. This stability criterion
asserts that a primary equilibrium state becomes unstable at a threshold where other nearby equilibrium states exist. In
the bifurcation theory, the existence of an incremental solution is just a necessary condition for the primary equilibrium
state at the threshold to be a bifurcation point. However, this condition is not suﬃcient. An extensive nonlinear bifurcation
analysis [60,74] is needed to verify the existence of the bifurcation solution branches, and therefore the adjacent equilibrium
states.
67
Yang & Sharma: A tutorial on the stability and bifurcation Page 68
In contrast to the implicit form, we consider the explicit forms of the bifurcation analysis of the
system shown in Fig. 16. Consider the equilibrium equations (273). The inﬁnitesimal increments of the
variables λ1, λ2, and ˜D are δλ1, δλ2, and δ ˜D, respectively. With the Taylor series expansion at the
equilibrium solutions (λ1, λ2, ˜D), we have
∂W
∂λ1
=
∂W
∂λ1 e
+
∂2
W
∂λ2
1 e
δλ1 +
∂2
W
∂λ1∂λ2 e
δλ2 +
∂2
W
∂λ1∂ ˜D e
δ ˜D + o(δλ1, δλ2, δ ˜D), (286a)
∂W
∂λ2
=
∂W
∂λ2 e
+
∂2
W
∂λ2∂λ1 e
δλ1 +
∂2
W
∂λ2
2 e
δλ2 +
∂2
W
∂λ2∂ ˜D e
δ ˜D + o(δλ1, δλ2, δ ˜D), (286b)
∂W
∂ ˜D
=
∂W
∂ ˜D e
+
∂2
W
∂ ˜D∂λ1 e
δλ1 +
∂2
W
∂ ˜D∂λ2 e
δλ2 +
∂2
W
∂ ˜D2
e
δ ˜D + o(δλ1, δλ2, δ ˜D). (286c)
Substituting Eq. (278) into the equilibrium equations (273) and ignoring the higher order terms, we have
∂W
∂λ1 e
+
∂2
W
∂λ2
1 e
δλ1 +
∂2
W
∂λ1∂λ2 e
δλ2 +
∂2
W
∂λ1∂ ˜D e
δ ˜D −
P1
L2L3
= 0, (287a)
∂W
∂λ2 e
+
∂2
W
∂λ2∂λ1 e
δλ1 +
∂2
W
∂λ2
2 e
δλ2 +
∂2
W
∂λ2∂ ˜D e
δ ˜D −
P2
L1L3
= 0, (287b)
∂W
∂ ˜D e
+
∂2
W
∂ ˜D∂λ1 e
δλ1 +
∂2
W
∂ ˜D∂λ2 e
δλ2 +
∂2
W
∂ ˜D2
e
δ ˜D −
Φ
L3
= 0. (287c)
Since the equilibrium solutions of λ1, λ2, and ˜D satisfy the equilibrium equations (273), that is
∂W
∂λ1 e
−
P1
L2L3
= 0,
∂W
∂λ2 e
−
P2
L1L3
= 0,
∂W
∂ ˜D e
−
Φ
L3
= 0, (288)
then Eq. (287) reduces to a set of simultaneous linear equations that can be expressed as a single matrix
equation, or, written out in full, as











∂2
W
∂λ2
1 e
∂2
W
∂λ1∂λ2 e
∂2
W
∂λ1∂ ˜D e
∂2
W
∂λ2∂λ1 e
∂2
W
∂λ2
2 e
∂2
W
∂λ2∂ ˜D e
∂2
W
∂ ˜D∂λ1 e
∂2
W
∂ ˜D∂λ2 e
∂2
W
∂ ˜D2
e





















δλ1
δλ2
δ ˜D










= 0. (289)
If the determinant of the 3×3 matrix in Eq. (289) is non-zero, the simultaneous linear equations only
have the trivial solution δλ1 = δλ2 = δ ˜D = 0. If the determinant is zero, there exists an inﬁnite number
of solutions of δλ1, δλ2, and δ ˜D. Thus, the necessary condition for the existence of incremental solutions
is a zero determinant of the 3 × 3 matrix in Eq. (289). Consider the free energy function (278) of an
ideal dielectric elastomer. By Eq. (289), the bifurcation condition of the electrostatic system shown in
68
Yang & Sharma: A tutorial on the stability and bifurcation Page 69
Fig. 16 at equilibrium (288) is given by
µ(1 + 3λ−4
1 λ−2
2 ) +
3 ˜D2
ε
λ−4
1 λ−2
2 2µλ−3
1 λ−3
2 +
2 ˜D2
ε
λ−3
1 λ−3
2 −
2 ˜D
ε
λ−3
1 λ−2
2
2µλ−3
1 λ−3
2 +
2 ˜D2
ε
λ−3
1 λ−3
2 µ(1 + 3λ−2
1 λ−4
2 ) +
3 ˜D2
ε
λ−2
1 λ−4
2 −
2 ˜D
ε
λ−2
1 λ−3
2
−
2 ˜D
ε
λ−3
1 λ−2
2 −
2 ˜D
ε
λ−2
1 λ−3
2
1
ε
λ−2
1 λ−2
2
= 0. (290)
5.2 Avoiding the pull-in instability of a dielectric elastomer ﬁlm
Inspired by the work by Zhao and Suo [38], Yang et al. [40] concocted a set of simple, experimentally
implementable, conditions that render the dielectric elastomer ﬁlm impervious to pull-in instability for
all practical loading conditions. Here we take the example in the work [40] to discuss electromechanical
instability and bifurcation.
5.2.1 Stability analysis by using the energy method
Figure 19: Schematic of a deformed ﬁlm of a dielectric elastomer. The ﬁlm can be extended/compressed
between two well-lubricated (with rollers), rigid plates by means of a controlled displacement (or the
stretch λ1) in the X1 direction. A dead load P2 is applied in the X2 direction and a voltage Φ is applied
in the thickness direction of the ﬁlm that is coated with two compliant electrodes. The controlled stretch
λ1, the dead load P2 and the voltage Φ deform the ﬁlm from L1, L2, and L3 to λ1L1, λ2L2, and λ3L3,
as well as induce an electric charge of magnitude Q on either electrode [40].
The free energy G of the system shown in Fig. 19 is [40]
G
L1L2L3
= W(λ1, λ2, ˜D) −
P2
L1L3
λ2 −
Φ
L3
˜D, (291)
where W(λ1, λ2, ˜D) is the free energy function. In contrast to three variables in Eq. (271), there are
only two variables λ2 and ˜D in Eq. (291). Here the reduction of variables from three to two is due to
69
Yang & Sharma: A tutorial on the stability and bifurcation Page 70
the displacement-controlled boundary condition introduced.
With the small variations δλ2 and δ ˜D, the variation of the free energy (291) is
∆G
L1L2L3
=
∂W
∂λ2
−
P2
L1L3
δλ2 +
∂W
∂ ˜D
−
Φ
L3
δ ˜D +
1
2
∂2
W
∂λ2
2
(δλ2)2
+
1
2
∂2
W
∂ ˜D2
(δ ˜D)2
+
∂2
W
∂λ2∂ ˜D
δλ2δ ˜D. (292)
In equilibrium, the ﬁrst variation of the free energy is zero, we have
∂W
∂λ2
−
P2
L1L3
= 0 and
∂W
∂ ˜D
−
Φ
L3
= 0. (293)
From the principle of minimum energy, the stability analysis is related to the Hessian matrix
H =





∂2
W
∂λ2
2
∂2
W
∂λ2∂ ˜D
∂2
W
∂ ˜D∂λ2
∂2
W
∂ ˜D2





. (294)
A stable equilibrium state requires that the Hessian matrix (294) must be positive deﬁnite, that is,
all the principal minors of H are positive in equilibrium. The 1 × 1 principal minors are
∂2
W
∂λ2
2
,
∂2
W
∂ ˜D2
, (295)
and the only 2 × 2 principal minor is
det H =
∂2
W
∂λ2
2
∂2
W
∂λ2∂ ˜D
∂2
W
∂ ˜D∂λ2
∂2
W
∂ ˜D2
. (296)
Again, if all the principal minors in Eqs. (295) and (296) are positive, the equilibrium state (293) is
stable. Consider the free energy function (278). The equilibrium equations (293) are
µ λ2 − 1 +
˜D2
εµ
λ−2
1 λ−3
2 =
P2
L1L3
and
˜D
ε
λ−2
1 λ−2
2 =
Φ
L3
. (297)
The Hessian matrix (294) is
H =







µ 1 + 3 1 +
˜D2
εµ
λ−2
1 λ−4
2 −
2 ˜D
ε
λ−2
1 λ−3
2
−
2 ˜D
ε
λ−2
1 λ−3
2
1
ε
λ−2
1 λ−2
2







. (298)
The 1 × 1 principle minors in Eq. (298) are the diagonal entries and they are always positive no
70
Yang & Sharma: A tutorial on the stability and bifurcation Page 71
matter the state is in equilibrium or not. The determinant of the Hessian in Eq. (298) is
det H =
µ
ε
λ−2
1 λ−3
2 4λ−2
1 λ−3
2 + λ2 − 1 +
˜D2
εµ
λ−2
1 λ−3
2 . (299)
By Eq. (297)1, we have
λ2 − 1 +
˜D2
εµ
λ−2
1 λ−3
2 =
1
µ
P2
L1L3
. (300)
Then we ﬁnd that the determinant at the equilibrium state is always positive
det H =
µ
ε
λ−4
1 λ−6
2 4 +
1
µ
P2
L1L3
λ2
1λ3
2 > 0 (301)
due to the positive stretches λ1 and λ2 and the positive applied dead load P2. Thus, a positive deﬁnite
Hessian matrix in this example indicates that the equilibrium state is stable.
5.2.2 Bifurcation analysis
The bifurcation analysis here is similar to that of the example in Sec. 5.1.2. Consider the equilibrium
equations (293). The inﬁnitesimal increments of the variables λ2 and ˜D are δλ2 and δ ˜D, respectively.
The Taylor series expansion gives
∂W
∂λ2
=
∂W
∂λ2 e
+
∂2
W
∂λ2
2 e
δλ2 +
∂2
W
∂λ2∂ ˜D e
δ ˜D + o(δλ2, δ ˜D), (302a)
∂W
∂ ˜D
=
∂W
∂ ˜D e
+
∂2
W
∂ ˜D∂λ2 e
δλ2 +
∂2
W
∂ ˜D2
e
δ ˜D + o(δλ2, δ ˜D). (302b)
Substituting Eq. (302) into the equilibrium equations (293) and ignoring the higher order terms, we have
∂W
∂λ2 e
+
∂2
W
∂λ2
2 e
δλ2 +
∂2
W
∂λ2∂ ˜D e
δ ˜D −
P2
L1L3
= 0, (303a)
∂W
∂ ˜D e
+
∂2
W
∂ ˜D∂λ2 e
δλ2 +
∂2
W
∂ ˜D2
e
δ ˜D −
Φ
L3
= 0. (303b)
Together with the equilibrium equations
∂W
∂λ2 e
−
P2
L1L3
= 0 and
∂W
∂ ˜D e
−
Φ
L3
= 0, (304)
we encounter a set of simultaneous linear equations that is expressed as a single matrix equation





∂2
W
∂λ2
2 e
∂2
W
∂λ2∂ ˜D e
∂2
W
∂ ˜D∂λ2 e
∂2
W
∂ ˜D2
e










δλ2
δ ˜D





= 0. (305)
71
Yang & Sharma: A tutorial on the stability and bifurcation Page 72
If the determinant of the 2 × 2 matrix in Eq. (305) is non-zero, namely
∂2
W
∂λ2
2 e
∂2
W
∂λ2∂ ˜D e
∂2
W
∂ ˜D∂λ2 e
∂2
W
∂ ˜D2
e
= 0, (306)
the simultaneous linear equations only have the trivial solution δλ2 = δ ˜D = 0. This means that there is
no bifurcation from the primary equilibrium state (solutions of Eq. (304)). Consider the example of an
ideal dielectric elastomer in Eqs. (297) and (298). The determinant of the Hessian matrix at equilibrium
is
µ 1 + 3 1 +
˜D2
εµ
λ−2
1 λ−4
2 −
2 ˜D
ε
λ−2
1 λ−3
2
−
2 ˜D
ε
λ−2
1 λ−3
2
1
ε
λ−2
1 λ−2
2
=
µ
ε
λ−4
1 λ−6
2 4 +
1
µ
P2
L1L3
λ2
1λ3
2 = 0, (307)
thus there is no bifurcation.
6 Stability of the homogeneous deformation of soft dielectrics:
an alternative energy formulation
In this section, we revisit the electromechanical instability problem presented in the work [38] but using
the energy formulation of nonlinear electroelasticity in terms of the deformation and the polarization
[62,63]. The procedure of the stability analysis can also be found in a recent paper by Alameh et al. [124]
which discussed a related problem of magneto-electro-mechanical instability.
6.1 A dielectric ﬁlm
Consider the homogeneous deformation of a thin ﬁlm of soft dielectrics subject to two pairs of equal inplane
forces (dead loads) and an electric voltage in the thickness direction. The dielectric ﬁlm deforms
from the original dimension LX × LY × LZ to λXLX × λY LY × λZLZ (see Fig. 20), such that
x = λXX, y = λY Y, z = λZZ, (308)
and the deformation gradient is
F := diag (λX, λY , λZ), (309)
where λX, λY , and λZ are constant stretches. Then, the Jacobian is
J = det F = λXλY λZ (310)
72
Yang & Sharma: A tutorial on the stability and bifurcation Page 73
and the inverses of F and FT
are
F−1
= F−T
:= diag (λ−1
X , λ−1
Y , λ−1
Z ). (311)
(a) (b)
Figure 20: Schematic of the deformation of a soft dielectric ﬁlm subjected to an electric voltage Φ and
two pairs of dead loads FX and FY . The dielectric ﬁlm is coated with two compliant electrodes on
the upper and bottom surfaces. (a) Undeformed with dimension LX × LY × LZ. (b) Deformed with
dimension λXLX × λY LY × λZLZ.
In the homogeneous deformation (308), the nominal polarization ˜P and the nominal electric ﬁeld ˜E
in Eq. (152) are uniform and only in the Z−direction, such that
˜P = − ˜PeZ, ˜E = − ξ = −
∂ξ
∂Z
eZ = −
Φ
LZ
eZ = − ˜E0eZ. (312)
With Eqs. (311) and (312), the true electric ﬁeld in Eq. (152) is
E = −F−T
ξ = −λ−1
Z
˜E0eZ, (313)
and then the nominal electric displacement in Eq. (156) reads
˜D = F−1
(ε0JE + ˜P) = −λ−1
Z (ε0λXλY
˜E0 + ˜P)eZ. (314)
The total free energy of the system shown in Fig. 20 is [63,125]
F(λX, λY , λZ, ˜P) = LXLY LZΨ(λX, λY , λZ, ˜P) − (FXλXLX + FY λY LY )
+ LXLY LZ
ε0
2
λ−1
Z λXλY
˜E2
0 − LXLY Φλ−1
Z (ε0λXλY
˜E0 + ˜P). (315)
Divided by the volume LXLY LZ, the free energy density reads
f(λX, λY , λZ, ˜P) = Ψ(λX, λY , λZ, ˜P) − SXλX − SY λY −
ε0
2
λ−1
Z λXλY
˜E2
0 − λ−1
Z
˜E0
˜P, (316)
where SX = FX/(LY LZ) and SY = FY /(LZLX). The internal energy function Ψ in Eq. (316) consists
73
Yang & Sharma: A tutorial on the stability and bifurcation Page 74
of two parts:
Ψ(λX, λY , λZ, ˜P) = We
(λX, λY , λZ) +
˜P2
2J(ε − ε0)
. (317)
The ﬁrst term We
on the RHS of Eq. (317) represents the elastic strain-energy function, whilst the
second term reﬂects the usual linear dielectric behavior, i.e., the permittivity ε of the dielectric elastomer
is independent of the deformation. Next, we discuss the strain-energy functions for both compressible
and incompressible materials.
6.1.1 Compressible soft dielectrics
The condition for the vanishing of the ﬁrst variation of Eq. (316) gives the equilibrium equations
∂f
∂λX
= 0,
∂f
∂λY
= 0,
∂f
∂λZ
= 0,
∂f
∂ ˜P
= 0. (318)
The Hessian matrix of Eq. (316) is
H =
















∂2
f
∂λ2
X
∂2
f
∂λX∂λY
∂2
f
∂λX∂λZ
∂2
f
∂λX∂ ˜P
∂2
f
∂λY ∂λX
∂2
f
∂λ2
Y
∂2
f
∂λY ∂λZ
∂2
f
∂λY ∂ ˜P
∂2
f
∂λZ∂λX
∂2
f
∂λZ∂λY
∂2
f
∂λ2
Z
∂2
f
∂λZ∂ ˜P
∂2
f
∂ ˜P∂λX
∂2
f
∂ ˜P∂λY
∂2
f
∂ ˜P∂λZ
∂2
f
∂ ˜P2
















. (319)
For isotropic compressible materials, the strain-energy function, for example, can be written as [63,
125,126]:
We
(λX, λY , λZ) =
µ
2
J−2/3
(λ2
X + λ2
Y + λ2
Z) − 3 +
κ
2
(J − 1)2
, (320)
where µ and κ are the shear and bulk moduli, respectively.
Substituting the free energy density (316), with Eqs. (317) and (320), into (318), we can obtain four
algebraic equations of four variables λX, λY , λZ, ˜P. Solutions of these four algebraic equations give
the equilibrium states (λX, λY , λZ, ˜P) of the dielectric ﬁlm shown in Fig. 20. Then the stability of the
equilibrium states can be veriﬁed by evaluating the positive deﬁniteness of the Hessian matrix (319) at
equilibrium.
6.1.2 Incompressible soft dielectrics
In contrast to compressible materials, the constraint of incompressibility requires that
J = det F = λXλY λZ = 1. (321)
i): Lagrange multiplier method
74
Yang & Sharma: A tutorial on the stability and bifurcation Page 75
The Lagrange multiplier approach is often used to address the extremum problem subjected to a
constraint, for example, the constraint of incompressibility [57, 81, 86]. Consider a Lagrange multiplier
La. The internal energy function Ψ in Eq. (317) is modiﬁed as
¯Ψ(λX, λY , λZ, ˜P, La) = Ψ(λX, λY , λZ, ˜P) − La(λXλY λZ − 1) (322)
and the free energy density f(λX, λY , λZ, ˜P) in Eq. (316) becomes
¯f(λX, λY , λZ, ˜P, La) = f(λX, λY , λZ, ˜P) − La(λXλY λZ − 1). (323)
Then, the governing equations are
∂ ¯f
∂λX
= 0,
∂ ¯f
∂λY
= 0,
∂ ¯f
∂λZ
= 0,
∂ ¯f
∂ ˜P
= 0,
∂ ¯f
∂La
= 0. (324)
ii): Reduced energy function method
The constraint of incompressibility (321) gives the following relation between the stretches
λZ =
1
λXλY
. (325)
Then, the internal energy function (317) can be deﬁned by a reduced form
ˆΨ(λX, λY , ˜P) = Ψ(λX, λY ,
1
λXλY
, ˜P) (326)
and the free energy density (316) becomes
ˆf(λX, λY , ˜P) = ˆΨ(λX, λY , ˜P) − SXλX − SY λY −
ε0
2
λ2
Xλ2
Y
˜E2
0 − λXλY
˜E0
˜P. (327)
The zero ﬁrst variation of Eq. (327) gives the equilibrium equations
∂ ˆf
∂λX
(λX, λY , ˜P) = 0,
∂ ˆf
∂λY
(λX, λY , ˜P) = 0,
∂ ˆf
∂ ˜P
(λX, λY , ˜P) = 0. (328)
The related Hessian matrix of Eq. (327) is
ˆH =











∂2 ˆf
∂λ2
X
∂2 ˆf
∂λX∂λY
∂2 ˆf
∂λX∂ ˜P
∂2 ˆf
∂λY ∂λX
∂2 ˆf
∂λ2
Y
∂2 ˆf
∂λY ∂ ˜P
∂2 ˆf
∂ ˜P∂λX
∂2 ˆf
∂ ˜P∂λY
∂2 ˆf
∂ ˜P2











. (329)
75
Yang & Sharma: A tutorial on the stability and bifurcation Page 76
6.1.3 Example: a ﬁlm of incompressible neo-Hookean dielectric
We take the incompressible neo-Hookean dielectric as an example to illustrate the above analysis. Thus
the internal energy function (317), by the relation λZ = 1/(λXλY ) and the reduced energy function
(326), may be rewritten as
ˆΨ(λX, λY , ˜P) =
µ
2
(λ2
X + λ2
Y +
1
λ2
Xλ2
Y
− 3) +
˜P2
2(ε − ε0)
. (330)
Then the reduced free energy density (327) becomes
ˆf(λX, λY , ˜P) =
µ
2
(λ2
X +λ2
Y +
1
λ2
Xλ2
Y
−3)+
˜P2
2(ε − ε0)
−SXλX −SY λY −
ε0
2
λ2
Xλ2
Y
˜E2
0 −λXλY
˜E0
˜P. (331)
By Eq. (331), the governing equations (328) are
µ(λX −
1
λ3
Xλ2
Y
) − SX − ε0λXλ2
Y
˜E2
0 − λY
˜E0
˜P = 0, (332a)
µ(λY −
1
λ2
Xλ3
Y
) − SY − ε0λ2
XλY
˜E2
0 − λX
˜E0
˜P = 0, (332b)
˜P
ε − ε0
− λXλY
˜E0 = 0, (332c)
and the Hessian matrix (329) reads
ˆH =












µ(1 +
3
λ4
Xλ2
Y
) − ε0λ2
Y
˜E2
0 µ
2
λ3
Xλ3
Y
− 2ε0λXλY
˜E2
0 − ˜E0
˜P −λY
˜E0
µ(1 +
3
λ2
Xλ4
Y
) − ε0λ2
X
˜E2
0 −λX
˜E0
sym
1
ε − ε0












. (333)
In the following, we will discuss the solutions of Eqs. (332a)−(332c), and then study the stability of
the dielectric ﬁlm by checking the positive deﬁniteness of the Hessian matrix (333) at equilibrium.
The equilibrium equation (332c) directly gives the relation
˜P = (ε − ε0)λXλY
˜E0. (334a)
By Eq. (334a), Eqs. (332a) and (332b) become
µ(λX −
1
λ3
Xλ2
Y
) − SX − ελXλ2
Y
˜E2
0 = 0, (334b)
µ(λY −
1
λ2
Xλ3
Y
) − SY − ελ2
XλY
˜E2
0 = 0, (334c)
76
Yang & Sharma: A tutorial on the stability and bifurcation Page 77
and the Hessian matrix (333) is recast as
ˆH =












µ(1 +
3
λ4
Xλ2
Y
) − ε0λ2
Y
˜E2
0 µ
2
λ3
Xλ3
Y
− (ε0 + ε)λXλY
˜E2
0 −λY
˜E0
µ(1 +
3
λ2
Xλ4
Y
) − ε0λ2
X
˜E2
0 −λX
˜E0
sym
1
ε − ε0












. (335)
Note that Eqs. (334b) and (334c) are the same as equations (6a) and (6b) in the work of Zhao and
Suo [38] with appropriate notation changes. For example, the subscripts ‘1’ and ‘2’ in the work [38]
correspond to the subscripts ‘X’ and ‘Y ’ in this tutorial. If we replace ˜D in (6a) and (6b) by ˜E in (6c)
in their work, then we get exactly the same two algebraic equations as Eqs. (334b) and (334c). This
equivalence serves as a direct validation of the alternative formulation presented in this tutorial. In the
following, we will use some special but easy cases to illustrate the behavior of a dielectric ﬁlm.
Consider the case of a uniaxial prestress, i.e., SY = 0. The equilibrium equation (334c) can be recast
as
λY = [λ2
X − λ4
X( ˜E0 ε/µ)2
]−1/4
. (336)
Substituting Eq. (336) into the equilibrium equation (334b), we get the following algebraic equation for
λX:
λX −
1 − λ2
X( ˜E0 ε/µ)2
λ2
X
−
( ˜E0 ε/µ)2
1 − λ2
X( ˜E0 ε/µ)2
=
SX
µ
. (337)
The equilibrium solution of λX in Eq. (337) for a given pair of (SX, ˜E0) is shown in Fig. 21(a).
Without the applied dead loads, i.e., SX = SY = 0, the nominal electric ﬁeld (prior to the onset of
instability) induces the ﬁlm to expand but not signiﬁcantly. When the nominal electric ﬁeld increases to
the threshold (i.e., the cross in each curve), pull-in instability occurs and the ﬁlm expands dramatically
without an appreciable increase in the electric ﬁeld. Finally, the thinned ﬁlm is damaged by electric
breakdown. In addition, the dead load SX expands the ﬁlm signiﬁcantly and the electric ﬁeld assists the
expansion, and pull-in instability occurs when the electric ﬁeld increases to the threshold. In contrast,
the threshold of the nominal electric ﬁeld ˜E0 decreases as the dead load SX increases.
To investigate electric actuation, we deﬁne the prestretch in the X−direction as λp
X, which is the
stretch due to the dead load SX in the absence of the voltage. Explicitly, λp
X is the solution of the
algebraic equation (337) with a zero electric ﬁeld ˜E0 = 0. In contrast to the prestretch λp
X, the actuation
stretch in the X−direction is deﬁned as λX/λp
X. Figure 21(b) shows the variation of the actuation
stretch λX/λp
X with respect to the nominal electric ﬁeld ˜E0 under several uniaxial dead loads SX. The
maximum actuation stretch in each curve is the critical stretch at the occurrence of pull-in instability
(marked by a cross). Without the dead load SX = 0, the ﬁlm can attain the largest maximum actuation
77
Yang & Sharma: A tutorial on the stability and bifurcation Page 78
stretch. By increasing the dead load SX, the maximum actuation stretch decreases. Note that even the
maximum actuation stretch decreases, the total stretch λX does increase as the increasing of the dead
load before the onset of pull-in instability. The diﬀerence between the total stretch λX and the actuation
stretch λX/λp
X can be observed by comparing Figs. 21(a) and 21(b).
(a) (b)
Figure 21: Deformation of the dielectric ﬁlm subject to an electric voltage and a uniaxial prestress
(SY = 0). (a) Stretch λX vs. the nominal electric ﬁeld under several forces SX/µ. (b) Actuation stretch
λX/λp
X vs. the nominal electric ﬁeld under several forces SX/µ. The critical point for the onset of pull-in
instability is denoted by a cross in each curve.
6.2 A dielectric disc
Another tractable example that we consider is that of a dielectric disc. To that end, we consider a
circular dielectric plate with radius Ra and thickness Ha in the undeformed state (see Fig. 22). Taking
the cylindrical coordinates (R, Θ, Z) with an orthonormal basis (eR, eΘ, eZ), the domain of the dielectric
disc in the reference conﬁguration is represented by
ΩR = (R, Θ, Z) ∈ R3
: 0 ≤ R ≤ Ra, 0 ≤ Θ < 2π, 0 ≤ Z ≤ Ha . (338)
6.2.1 Large deformation
Subject to both an applied voltage Φ in the thickness direction (the Z−direction) and a surrounding
uniform traction force FR, the dielectric disc expands its in-plane area and decreases its thickness from
Ha to λZHa. Here the stretch λZ in the thickness direction is a constant that is independent of the
position. We assume that the dielectric disc admits the cylindrical deformation
r = r(R), θ = Θ, z = λZZ. (339)
The deformation gradient in the cylindrical coordinates is
F := diag (
dr
dR
,
r
R
, λZ). (340)
78
Yang & Sharma: A tutorial on the stability and bifurcation Page 79
(a) (b)
Figure 22: Schematic of the deformation of a dielectric disc subjected to an electric voltage Φ and a
surrounding dead load FR. The disc is coated with two compliant electrodes on the upper and bottom
surfaces. (a) Undeformed disc with radius Ra and thickness Ha. (b) Deformed disc with thickness λZHa.
Incompressibility requires that
J = det F =
dr
dR
r
R
λZ = 1. (341)
Integrating Eq. (341) with respect to R and with r(0) = 0, we have
r(R) =
R
√
λZ
. (342)
Thus, the deformation gradient (340) becomes
F := diag (
1
√
λZ
,
1
√
λZ
, λZ). (343)
The thinning of the disc shown in Fig. 22 is similar to the deformation in Fig. 20. By the cylindrical
deformation (339), the equations of the polarization, the electric ﬁeld, and the electric displacement are
the same as Eqs. (312)−(314).
With the energy formulation in the work [63,125] and the deformation gradient (343), the free energy
of the dielectric disc in Fig. 22 is written as
F(λZ, ˜P) = πR2
aHaΨ(λZ, ˜P) −
FRRa
√
λZ
+ πR2
aHa
ε0
2
λ−2
Z
˜E2
0 − πR2
aΦλ−1
Z (ε0λ−1
Z
˜E0 + ˜P). (344)
Divided by the volume πR2
aHa, the free energy density reads
f(λZ, ˜P) = Ψ(λZ, ˜P) −
SR
√
λZ
−
ε0
2
λ−2
Z
˜E2
0 − λ−1
Z
˜E0
˜P, (345)
where SR = FR/(πRaHa). Consider incompressible neo-Hookean dielectrics. By Eq. (343), the internal
energy function reads
Ψ(λZ, ˜P) =
µ
2
(λ2
Z + 2λ−1
Z − 3) +
˜P2
2(ε − ε0)
(346)
and the free energy density (345) becomes
f(λZ, ˜P) =
µ
2
(λ2
Z + 2λ−1
Z − 3) +
˜P2
2(ε − ε0)
−
SR
√
λZ
−
ε0
2
λ−2
Z
˜E2
0 − λ−1
Z
˜E0
˜P. (347)
79
Yang & Sharma: A tutorial on the stability and bifurcation Page 80
The zero ﬁrst variation of Eq. (347) gives
µ(λZ − λ−2
Z ) +
1
2
SRλ
−3/2
Z + ε0λ−3
Z
˜E2
0 + λ−2
Z
˜E0
˜P = 0, (348a)
˜P
ε − ε0
− λ−1
Z
˜E0 = 0. (348b)
Equation (348b) directly gives the relation
˜P = (ε − ε0)λ−1
Z
˜E0. (349)
Substituting Eq. (349) into Eq. (348a), we have
µ(λZ − λ−2
Z ) +
1
2
SRλ
−3/2
Z + ελ−3
Z
˜E2
0 = 0. (350)
Equation (350) is an algebraic equation of λZ with two controlled parameters SR and ˜E0. For a given
pair of (SR, ˜E0), there may exist many cases of the roots of λZ in Eq. (350), for example, multiple real
and complex roots. The detailed discussion of the properties of roots of λZ in Eq. (350) is not addressed
here. The solution for realistic physical situations of the stretch λZ, a real value and 0 < λZ ≤ 1, for a
given pair of (SR, ˜E0) is shown in Fig. 23.
Figure 23(a) shows the equilibrium solution of the stretch λZ in Eq. (350) for a given pair of (SR, ˜E0).
Without the dead load SR = 0, the stretch λZ decreases from 1 as the nominal electric ﬁeld ˜E0 increases.
With the increase of the electric ﬁeld, the disc decreases its thickness to the threshold at which pull-in
instability occurs. The onset of pull-in instability is at the peak of each curve. Later we will show that
the peak corresponds to a zero determinant of the Hessian matrix at equilibrium. With the increase
of the dead load SR, the peak in each curve decreases, which means that it is easier to induce pull-in
instability at a larger dead load. After the peak, the disc expands and the thickness decreases rapidly
and simultaneously until the onset of electric breakdown.
In contrast to the deformation in the thickness direction, Fig. 23(b) shows the deformation in the
radial direction with the stretch 1/
√
λZ. The curves in Figs. 23(b) and 21(a) have the same trends as
the increase of the applied dead load and the electric ﬁeld.
Figure 23(c) shows the variation of the stretch λZ in the thickness direction with respect to the dead
load SR for varying electric ﬁelds. Interestingly, the curve with a zero electric ﬁeld has no peak, and the
applied dead load increases monotonically with the decrease of the thickness (λZ).
Figure 23(d) shows the variation of the stretch 1/
√
λZ of the radial deformation with respect to the
dead load SR under several electric ﬁelds. At a zero dead load SR = 0, it is obvious that the stretches
correspond to diﬀerent electric ﬁelds are almost the same and are close to 1, which indicates that the
80
Yang & Sharma: A tutorial on the stability and bifurcation Page 81
(a) (b)
(c) (d)
Figure 23: Deformation of the dielectric disc subject to the electric voltage and the surrounding force.
(a) Stretch λZ in the thickness direction vs. the nominal electric ﬁeld under several surrounding forces
SR/µ. (b) Stretch 1/
√
λZ in the radial direction vs. the nominal electric ﬁeld under several surrounding
forces SR/µ. (c) Stretch λZ in the thickness direction vs. the surrounding force SR/µ under several
nominal electric ﬁelds. (d) Stretch 1/
√
λZ in the radial direction vs. the surrounding force SR/µ under
several nominal electric ﬁelds. The critical point for the onset of pull-in instability is denoted by a cross
in each curve.
disc deformation is less insensitive to the electric ﬁeld compared to the dead load. However, the trend of
each curve changes signiﬁcantly with the increase of the applied dead load, especially the peak of each
curve. This implies that for a given electric ﬁeld, the increase of the applied dead load will eventually
make pull-in instability occur. For a larger electric ﬁeld, the peak of the curve occurs at a smaller dead
load.
All the curves in Fig. 23 simply emerge from the equilibrium equation (350), and a point on the curve
corresponds to an equilibrium state. The stability of the equilibrium states has not veriﬁed so far. When
we described the trend of each curve in Fig. 23, we claimed that pull-in instability occurs at the peak
of each curve and the disc is stable before the point reaches the peak. We will verify this claim in the
following stability analysis.
81
Yang & Sharma: A tutorial on the stability and bifurcation Page 82
6.2.2 Stability analysis and actuation strain
In this section, we discuss the stability and (pull-in) instability regions of the dielectric disc subjected to
both the dead load and the applied voltage. We also discuss the critical stretch and the critical actuation
stretch (the maximum actuation stretch) at which pull-in instability occurs.
The stability of the equilibrium solutions in Eqs. (349) and (350), see Fig. 23, can be determined by
examining the positive deﬁniteness of the Hessian matrix at equilibrium. By Eq. (349) and the second
variation of the free energy (347), we have the Hessian matrix
H =




µ(1 + 2λ−3
Z ) −
3
4
SRλ
−5/2
Z − (ε0 + 2ε)λ−4
Z
˜E2
0 λ−2
Z
˜E0
λ−2
Z
˜E0
1
ε − ε0



 (351)
with the determinant
det H =
1
ε − ε0
µ(1 + 2λ−3
Z ) −
3
4
SRλ
−5/2
Z − 3ελ−4
Z
˜E2
0 . (352)
The stability of the homogeneous deformation requires that the 2 × 2 Hessian matrix (351) should be
positive deﬁnite at the equilibrium solution (350), i.e., all the principal minors of Eq. (351) are positive
subject to the equilibrium equation (350). The 1 × 1 principal minors are the diagonal entries H11 and
H22, and the only 2 × 2 principal minor is the determinant (352). Since ε is greater than ε0, the entry
H22 is always positive. More importantly, compared to the entry H11, the determinant (352) is more
likely to become zero from positive for an equilibrium solution of λZ in Eq. (350). Therefore, in the case
of the homogeneous deformation of a dielectric disc, from the Hessian method the equilibrium solution
of λZ from Eq. (350) becomes unstable when the determinant (352) decreases to zero at the equilibrium
solution λZ, that is det H = 0, namely
µ(1 + 2λ−3
Z ) −
3
4
SRλ
−5/2
Z − 3ελ−4
Z
˜E2
0 = 0, (instability condition) (353)
where λZ is the solution of the equilibrium equation (350) for a given pair of (SR, ˜E0).
On each equilibrium curve in Fig. 23, the occurrence of pull-in instability is marked by a cross at
which both the equilibrium equation (350) and the instability condition (353) are satisﬁed. Through
the peak is determined by the solution of both Eqs. (350) and (353), it can be understood from the
conventional extremization problem of a scalar function with a scalar variable.
Again, Eqs. (350) and (353) together deﬁne the threshold (either the critical electric ﬁeld ˜Ec
0 or the
critical dead load) of pull-in instability. For example, for a given dead load SR, the solution of the
equilibrium equation (350) and the instability condition (353) is λZ and ˜Ec
0. For any electric ﬁeld with
a magnitude smaller than ˜Ec
0, there is no pull-in instability.
82
Yang & Sharma: A tutorial on the stability and bifurcation Page 83
(a) (b)
Figure 24: The stability and instability regions in the dead load - the electric ﬁeld plane. (a) The dead
load SR vs. the nominal electric ﬁeld ˜E0. (b) The dead load SR vs. the true electric ﬁeld E0. The critical
true electric ﬁeld of pull-in instability is essential for the analysis of electric breakdown of the dielectric
disc.
Figure 24 plots the stability and instability regions in the dead load - the electric ﬁeld plane. The
stability and instability regions are separated by a blue curve, which is the solution ˜Ec
0 of the equilibrium
equation (350) and the instability condition (353) for an applied dead load SR. For a given dead load
SR, if the applied electric ﬁeld is less than ˜Ec
0, there is no pull-in instability and the dielectric disc is
stable. When the electric ﬁeld increases to ˜Ec
0 at which the instability condition (353) is satisﬁed, pull-in
instability occurs. Then the dielectric disc expands its area rapidly, and ﬁnally, electric breakdown occurs
at a relatively high true electric ﬁeld.
(a) (b)
Figure 25: Eﬀects of the surrounding force SR/µ on (a) the prestretch λp
Z, the critical stretch λc
Z and
the actuation stretch λc
Z/λp
Z in the thickness direction and (b) the prestretch 1/ λp
Z, the critical stretch
1/ λc
Z and the actuation stretch λp
Z/ λc
Z in the radial direction.
The actuation stretch problem at equilibrium has been discussed in the dielectric ﬁlm problem in
83
Yang & Sharma: A tutorial on the stability and bifurcation Page 84
Fig. 21(b). A natural question to ask here is how large can the actuation stretch be achieved before the
occurrence of pull-in instability. In Fig. 21(b), we know that the maximum actuation stretch is obtained
at the peak point (the cross) of each curve. In this dielectric disc problem, we plot the variation of the
maximum actuation stretch with the increase of the dead load in Fig. 25. Also, the prestretch and the
total stretch are plotted for a direct comparison.
In Fig. 25(a), we plot the maximum actuation stretch λc
Z/λp
Z in the thickness direction (Z−direction).
Since both the dead load and the electric ﬁeld decrease the thickness of the disc, the prestretch λp
Z, the
critical stretch λc
Z, and the maximum actuation stretch λc
Z/λp
Z are all less than 1. With the increase
of the dead load, both the prestretch λp
Z and the critical stretch λc
Z decrease rapidly to zero. However,
the maximum actuation stretch λc
Z/λp
Z decreases slowly and then reaches a plateau as the dead load
increases further.
In contrast, we plot the maximum actuation stretch λp
Z/ λc
Z in the radial direction in Fig. 25(b).
Both the prestretch 1/ λp
Z and the critical stretch 1/ λc
Z increase monotonically with the increase of
the dead load. However, the maximum actuation stretch is almost constant around the value 1.5. This
suggests a useful guideline for the design of a dielectric disc actuator i.e., the maximum actuation stretch
is insensitive to the applied dead load.
6.2.3 Pull-in instability vs. electric breakdown
Figure 26: The competition between pull-in instability and electric breakdown in the dead load - the
true electric ﬁeld plane.
Electric breakdown occurs in the dielectric elastomer when the true electric ﬁeld reaches a relatively
high value beyond some threshold [5, 13, 40]. Consider a true electric ﬁeld EEB = 3 × 108
V m−1
at
84
Yang & Sharma: A tutorial on the stability and bifurcation Page 85
which electric breakdown (EB) occurs. Other material parameters used in the numerical calculations
are µ = 106
N m−2
and ε = 3.54 × 10−11
F m−1
. Then the dimensionless true electric ﬁeld is given by
EEB ε/µ = 1.785. Thus, the stability and (pull-in) instability regions in Fig. 24 are altered.
Since EEB ε/µ = 1.785 is a constant,14
the threshold of electric breakdown is just a horizontal line,
i.e., the true electric ﬁeld E0 ε/µ = 1.785 in Fig. 26. Above the horizontal line, electric breakdown
would occur in the dielectric disc. In addition, below the horizontal line, there is no electric breakdown,
and the disc can be either stable or vulnerable to pull-in instability. Thus, there exists a competition
between electric breakdown and pull-in instability. This is graphically depicted in Fig. 26.
7 Electromechanical cavitation and the snap-through instability
in soft hollow dielectric spheres
In this section, we study the nonlinear electromechanical behavior of soft hollow dielectric spheres. For
the purely mechanical case, the problem can be reduced to the well-studied phenomenon of cavitation in
nonlinear elasticity studied by Ball [127]. Another related problem is the inﬂation of a soft balloon within
the framework of nonlinear elasticity. The inﬂation process accompanies snap-through instability that
can be harnessed to achieve giant voltage-triggered deformation of dielectric balloons [128,129]. Nonlinear
coupling problem in spherical shells of dielectrics can also be found in other works [107, 130–134], for
example, the snap-through instability of a thick-wall electro-active balloon [135] and multilayer electroactive
spherical balloons [136].
7.1 Formulation
Consider a spherical dielectric shell with thickness H and inner radius Ra and outer radius Rb. Spherical
coordinates (R, Θ, Φ) with corresponding local base unit vectors (eR, eΘ, eΦ) will be employed. The
domain occupied by the undeformed spherical dielectric (see Fig. 27) is represented by
ΩR = (R, Θ, Φ) ∈ R3
: Ra ≤ R ≤ Rb, 0 ≤ Θ < π, 0 ≤ Φ < 2π . (354)
The boundary of ΩR is denoted by ∂ΩR consisting of the inner and outer surfaces:
∂ΩRa = {X ∈ ΩR : |X| = Ra} and ∂ΩRb
= {X ∈ ΩR : |X| = Rb} . (355)
14The true electric ﬁeld for the onset of electric breakdown may slightly vary with the deformation of a dielectric
elastomer. The constant true electric ﬁeld EEB is just an approximation used here for the illustration of the competition
between pull-in instability and electric breakdown.
85
Yang & Sharma: A tutorial on the stability and bifurcation Page 86
(a)
(b) (c)
Figure 27: Schematic of the deformation of a hollow sphere of soft dielectrics subjected to the dead loads
Sa and Sb on the inner and outer surfaces, and an applied voltage Φ. The inner and outer surfaces are
coated with compliant electrodes. (a) 3D view of an undeformed hollow sphere with inner radius Ra and
outer radius Rb. (b) 2D view of an undeformed hollow sphere. (c) 2D view of a deformed hollow sphere
with inner radius ra and outer radius rb.
7.1.1 Energy of the system
The free energy of the system consists of two parts [35,63]:
F(χ, ˜P) = Umech
(χ) + Eelct
(χ, ˜P). (356)
The term Umech
(χ) is the mechanical part
Umech
(χ) =
ΩR
We
( χ) −
∂ΩR
tR · χ (357)
while the term Eelct
(χ, ˜P) is the electric part, for linear dielectrics, it can be expressed as
Eelct
(χ, ˜P) =
ε0
2 ΩR
J F−T
ξ
2
+
∂ΩR
ξ − ε0JC−1
ξ + F−1 ˜P · N +
ΩR
|˜P|2
2J(ε − ε0)
, (358)
86
Yang & Sharma: A tutorial on the stability and bifurcation Page 87
where N is the outer unit normal to ∂ΩR. The Maxwell equation is
Div − 0JC−1
ξ + F−1 ˜P = 0 in ΩR (359)
and the electric boundary conditions are
ξ = ξa on ∂ΩRa
and ξ = ξb on ∂ΩRb
. (360)
7.1.2 First variation of the energy functional
Since the deformation gradient F = χ and the polarization ˜P are independent, the ﬁrst variation of
the free energy (356) with respect to ˜P and F, respectively, gives (see the Appendix in [35] for details)
˜P = −(ε − ε0)JF−T
ξ, Div T = 0 in ΩR, (361)
TN = tR on ∂ΩR, (362)
where the total nominal stress T is
T =
∂We
∂F
+ TM
− LaF−T
. (363)
Here La serves as the Lagrange multiplier, TM
denotes the Maxwell stress
TM
=
1
εJ
(F ˜D) ⊗ ˜D −
1
2εJ
|F ˜D|2
F−T
, (364)
and
˜D = −εJC−1
ξ. (365)
Note that the total nominal stress T and the Maxwell stress TM
can also be found in Eqs. (195a) and
(192), respectively. The polarization ˜P and the electric displacement ˜D are consistent with Eqs. (158)
and (157) for linear dielectrics. Now we have formed a boundary-value problem (359) − (362).
7.1.3 Spherical deformation
A deformation x = χ(X) is spherical if it has the following component form
r = r(R), ϑ = Θ, ϕ = Φ, (366)
where r, ϑ and ϕ are the spherical coordinates. The spherical deformation implies that the hollow sphere
also admits a spherical shape after the deformation (see Fig. 27).
In the spherical coordinates, the deformation gradient of the spherical deformation (366) has the
87
Yang & Sharma: A tutorial on the stability and bifurcation Page 88
following component form
F := diag (λ1, λ2, λ3), (367)
where
λ1 = r , λ2 = λ3 =
r
R
. (368)
In this section, a prime denotes the derivative with respect to R. Consequently,
F−1
= F−T
:= diag (λ−1
1 , λ−1
2 , λ−1
3 ), C−1
:= diag (λ−2
1 , λ−2
2 , λ−2
3 ). (369)
In the spherical deformation of a hollow dielectric sphere with an applied voltage on the inner and
outer surfaces, the electric potential ξ is only a function of R, i.e., ξ = ξ(R). And the boundary conditions
(360) can be expressed as ξa = ξ(Ra) and ξb = ξ(Rb). Hence, the gradient of the potential in the spherical
coordinates in the reference conﬁguration reads
ξ(R) = ξ eR. (370)
With Eqs. (369) and (370), the nominal electric displacement ˜D in Eq. (365), the Maxwell equation
(359), and the Maxwell stress TM
in Eq. (364) become
˜D = −εJλ−2
1 ξ eR, (371)
εJ (λ−2
1 ξ ) +
2λ−2
1 ξ
R
= 0, (372)
and
TM
:=
εJλ−2
1 (ξ )2
2
diag (λ−1
1 , −λ−1
2 , −λ−1
3 ). (373)
Invoking isotropy, the strain-energy function We
and the tensor
∂We
∂F
depend on F = χ through
the principal stretches λ1, λ2 and λ3 (see Eqs. (181) and (182)). The total nominal stress (363) is
represented by
T := diag (T11, T22, T33), (374)
where
Tii =
∂We
∂λi
− Laλ−1
i + TM
ii , i = 1, 2, 3 (no sum). (375)
By Eq. (374), the equilibrium equation (361)2 gives the only non-trivial equation
∂T11
∂R
+
2T11
R
−
T22 + T33
R
= 0, (376)
and the natural boundary conditions (362) read
T11 = Sa at R = Ra and T11 = Sb at R = Rb. (377)
88
Yang & Sharma: A tutorial on the stability and bifurcation Page 89
Note that λ2 = λ3 in Eq. (367) and then T22 = T33 in Eq. (376).
7.1.4 Ideal dielectric elastomer
By the strain-energy function of incompressible neo-Hookean solids in Eq. (184), Eq. (375) becomes
Tii = µλi − Laλ−1
i + TM
ii , i = 1, 2, 3 (no sum), (378)
and the unit Jacobian J = det F = 1 yields
r r2
R2
= λ1λ2λ3 = 1. (379)
Let us introduce a new variable
λ = λ2 = λ3 =
r(R)
R
. (380)
It follows from Eqs. (368) and (379) that
r = λ1 = λ−1
2 λ−1
3 = λ−2
. (381)
Then, the total nominal stress (378) is reduced to
T11 = µ +
ε(ξ )2
2
λ8
λ−2
− Laλ2
, T22 = T33 = µ −
ε(ξ )2
2
λ2
λ − Laλ−1
. (382)
The traction boundary conditions (377) become
µ +
ε(ξ )2
2
λ8
λ−2
− Laλ2
= Sa at R = Ra, (383)
µ +
ε(ξ )2
2
λ8
λ−2
− Laλ2
= Sb at R = Rb. (384)
By the chain-rule, we have
∂
∂R
=
∂λ
∂R
∂
∂λ
=
1
R
(λ−2
− λ)
∂
∂λ
. (385)
By Eqs. (385), (372) and (382), the equilibrium equation (376) can be reduced to
2µ(1 − λ−3
) − λ2 ∂La
∂λ
= 0. (386)
7.1.5 Solution of the boundary-value problem
We now have a boundary-value problem consisting of the equilibrium equation (386), the traction boundary
conditions (383) and (384), and the electric boundary conditions (360). Integration of Eq. (386) with
respect to λ gives
La = 2µ
1
4
λ−4
− λ−1
+ Ccons, (387)
89
Yang & Sharma: A tutorial on the stability and bifurcation Page 90
where Ccons is an integral constant. By Eq. (379), the integration of Eq. (372) gives
ξ = C0
r
r2
, ξ(R) = −
C0
r(R)
+ C1, (388)
where C0 and C1 are the integration constants. With the electric boundary conditions ξa = ξ(Ra) and
ξb = ξ(Rb), we have
C0 =
rarb
rb − ra
(ξb − ξa) =
rarb
rb − ra
Φ, C1 =
rbξb − raξa
rb − ra
. (389)
With Eqs. (379) and (381), ξ in Eq. (388)1 can be recast as ξ = C0λ−4
/R2
. By Eq. (387), the
boundary conditions (383) and (384) are equal to
1
2
µλ−4
a + 2µλ−1
a − Ccons +
εC2
0
2R4
a
λ−4
a = Saλ−2
a , (390)
1
2
µλ−4
b + 2µλ−1
b − Ccons +
εC2
0
2R4
b
λ−4
b = Sbλ−2
b . (391)
Deleting the constant Ccons in Eqs. (390) and (391), we ﬁnally have the relation
1
2
µ(λ−4
b − λ−4
a ) + 2µ(λ−1
b − λ−1
a ) +
εC2
0
2
λ−4
b
R4
b
−
λ−4
a
R4
a
= Sbλ−2
b − Saλ−2
a , (392)
where
C0 =
rarb
rb − ra
Φ =
λaλbRaRb
λbRb − λaRa
Φ, λa =
ra
Ra
and λb =
rb
Rb
=
3
r3
a + R3
b − R3
a
Rb
. (393)
7.2 Electromechanical cavitation and the snap-through instability
7.2.1 Ball’s classic cavitation problem
In Ball’s cavitation problem [127], there is no electric ﬁeld within the hollow sphere, that is C0 = 0 in
Eq. (392). Also, the traction force on the inner surface is zero, Sa = 0. Thus, Eq. (392) reduces to
Sb
2µ
= λb +
1
4
λ−2
b − λ−1
a +
1
4
λ−4
a λ2
b. (394)
With λa and λb in Eq. (393), Eq. (394) is given by
Sb
2µ
= ¯r3
a + 1 − ¯R3
a
1/3
+
1
4
¯r3
a + 1 − ¯R3
a
−2/3
− ¯r−1
a
¯Ra +
1
4
¯r−4
a
¯R4
a ¯r3
a + 1 − ¯R3
a
2/3
, (395)
where the nondimensional variables are
¯ra =
ra
Rb
, ¯Ra =
Ra
Rb
. (396)
In the following, some limiting cases based on the geometry and the material properties of the hollow
sphere are discussed. In the case of a suﬃciently small hollow sphere, i.e., ¯Ra → 0, at ﬁxed ¯ra, the third
90
Yang & Sharma: A tutorial on the stability and bifurcation Page 91
term on the RHS of Eq. (395) vanishes and the other two terms give
Sb
2µ
→
5/4 + ¯r3
a
(1 + ¯r3
a)2/3
, (397)
which implies Sb
2µ → 5
4 as ¯ra → 0. The limiting case can be found in the numerical plots in Fig. 29(a).
In Ball’s paper [127],
(a) (b)
Figure 28: Variation of the radius ra of the cavity with the applied stress Sb on the outer surface and the
applied electric ﬁeld E0 in the radial direction of the hollow sphere. The ﬁgures have been drawn for the
ratio Ra/Rb = 0.3. (a) Variation of the cavity radius with respect to the applied dead load under several
electric ﬁelds. (b) Variation of the cavity radius with the applied electric ﬁeld under several applied dead
loads. A snap-through instability phenomenon can be observed by increasing the applied electric ﬁeld
to the threshold.
(a) (b)
Figure 29: Variation of the radius ra of the cavity with an applied stress Sb on the outer surface of the
hollow sphere and an applied electric ﬁeld E0 in the radial direction. (a) Variation of the cavity radius
with the applied dead load under several ratios Ra/Rb without electric ﬁelds. (b) Variation of the cavity
radius with the applied dead load under several electric ﬁelds at a ratio Ra/Rb = 0.001.
91
Yang & Sharma: A tutorial on the stability and bifurcation Page 92
7.2.2 Eﬀect of the electric ﬁeld on cavitation and snap-through instability
Consider only the applied dead load on the outer surface and the applied voltage, i.e., Sa = 0 in Eq. (392).
Then, Eq. (392) reduces to
Sb
2µ
= λb +
1
4
λ−2
b − λ−1
a +
1
4
λ−4
a λ2
b +
εC2
0
4µ
λ−2
b
R4
b
−
λ−4
a λ2
b
R4
a
. (398)
By ¯Ra in Eq. (396) and the nominal electric ﬁeld ˜E0 = Φ/(Rb − Ra) = Φ/H, Eq. (398) becomes
Sb
2µ
= λb +
1
4
λ−2
b − λ−1
a +
1
4
λ−4
a λ2
b +
λ2
a
¯R2
a(1 − ¯Ra)2
(1 − λ−4
a
¯R−4
a λ4
b)
4(λb − λa
¯Ra)2
˜E0
µ/ε
2
. (399)
Thus, for a given dead load Sb/2µ, an applied electric ﬁeld ˜E0 ε/µ, and a given geometry of the
hollow sphere ¯Ra = Ra/Rb, Eq. (399) determines the deformed inner radius ¯ra = ra/Rb. The snapthrough
instability phenomenon can be observed by evaluating Eq. (399) numerically (see Fig. 28(b)).
A similar form of Eq. (399) can also be found in other works [135,136].
8 Post-buckling analysis: “Lyapunov-Schmidt-Koiter” approach
For most of the stability problems in nonlinear elasticity, the analytical solution of the stability condition
is extremely diﬃcult. Furthermore, even the existence of a solution is not assured. Beginning with the
pioneering work by Koiter [137], the asymptotic expansion technique that follows the Lyapunov-Schmidt
decomposition has become a powerful tool to analyze the bifurcation and the (initial) post-bifurcation of
perfect or imperfect systems [55,61,65,67,138–143]. This procedure of the analysis of stability is known
as the “Lyapunov-Schmidt-Koiter” approach.
In this tutorial, we will go through the basic idea of the “Lyapunov-Schmidt-Koiter” approach and
then incorporate it into the post-buckling analysis of a dielectric elastomer by simple examples. Though
the “Lyapunov-Schmidt-Koiter” approach has been fruitfully applied to study the initial post-bifurcation
stage and investigate purely mechanical instability, only a few works [144] have reported the use of the
LSK approach in a post-bifurcated solution of the electromechanical coupling problem.
We consider the case of only one generalized state variable u of a 1D system Ω1
= {x ∈ R : 0 ≤ x ≤ l},
i.e., u ∈ C2
([0, l]; Rn
). For example, u can be the displacement of deformable solid.15
Thus, the potential
energy P of the system is a functional of u : [0, l] → Rn
. In addition, P depends on the load parameter
applied to the system. Here we assume the load parameter is a single scalar variable λ ∈ R that
corresponds to the work done by the external load. The potential energy of the conservative system can
15In electrostatics of deformable media, there are always two state variables: one corresponds to the deformation and
the other corresponds to the electric quantity (the electric displacement [95], the electric ﬁeld [92] or the polarization [63]).
However, the central idea of the LSK approach for multiple state variables is the same. Accordingly, for simplicity, we
illustrate with just one generalized state variable u.
92
Yang & Sharma: A tutorial on the stability and bifurcation Page 93
be written as
P[u; λ]. (400)
At a given load parameter λ ∈ R, an equilibrium state u : [0, l] → Rn
requires that
δP[u; λ] =
∂P[u + τδu; λ]
∂τ τ=0
≡ P [u; λ]δu = 0 (401)
for all kinematically admissible variations δu ∈ C2
([0, l]; Rn
).
For convenience, we introduce the following notations of the derivatives of functionals [67]:
P [u; λ]u1 ≡
∂P[u + τ1u1; λ]
∂τ1 τ1=0
,
P [u; λ]u1u1 ≡ P [u; λ]u2
1 ≡
∂2
P[u + τ1u1; λ]
∂2τ1 τ1=0
,
P [u; λ]u1u2 ≡ P [u; λ]u2u1 ≡
∂2
P[u + τ1u1 + τ2u2; λ]
∂τ1∂τ2 τ1=τ2=0
,
P [u; λ]u1u2u3 ≡
∂3
P[u + τ1u1 + τ2u2 + τ3u3; λ]
∂τ1∂τ2∂τ3 τ1=τ2=τ3=0
,
...
P(n)
[u; λ]u1u2 · · · un ≡
∂(n)
P u +
n
i=1
τiui; λ
∂τ1∂τ2 · · · ∂τn
τ1=τ2=···=τn=0
,



(402)
where u ∈ C2
([0, l]; Rn
), u1, u2, u3, . . . , un ∈ C2
([0, l]; Rn
), and τ1, τ2, . . . , τn, λ ∈ R.
8.1 Bifurcation analysis
We assume that there exists a principal (trivial) branch u0 that varies smoothly with the load parameter
λ as the load increases from zero. The equilibrium of the principal branch u0, from Eq. (401), reads
δP[u0(λ); λ] ≡ P [u0(λ); λ]δu = 0. (403)
As the load parameter λ increases from zero to the threshold λc, there exists a bifurcated branch
u(λ) that bifurcates from the principal branch u0(λ), such that
u(λ) = u0(λ) + v(λ). (404)
The bifurcated branch u(λ) intersects the trivial branch u0(λ) at λc, namely
lim
λ→λc
v(λ) = 0. (405)
93
Yang & Sharma: A tutorial on the stability and bifurcation Page 94
The bifurcation buckling mode can be deﬁned by [67]
u1 = lim
λ→λc
v(λ)
v
, (406)
where · represents a suitable norm, for example, the inner product form u1 = u1, u1
1/2
. The norm
of the buckling mode u1 deﬁned in Eq. (406) is unity, i.e., u1 = 1. In this tutorial, we only consider
the case of a single buckling mode at λ = λc. A classic example of the single mode at the threshold
is the buckling of Euler’s column. In contrast, surface instability of a compressed elastic half-space is
the case of the multiple buckling modes [16,70,145], and the LSK approach is used to study the initial
post-bifurcation behavior of surface wrinkling of a compressed half-space [145] and a neo-Hookean bilayer
[146], etc.
Similar to Eq. (403), the equilibrium equation (404) of the bifurcated branch is
P [u0(λ) + v(λ); λ]δu = 0. (407)
The Taylor series expansion of Eq. (407) is
P [u0; λ]δu + P [u0; λ]vδu +
1
2
P [u0; λ]v2
δu + · · · = 0. (408)
The ﬁrst term vanishes from Eq. (403). In the neighborhood of λc, |λ − λc| 1, dividing the
remaining terms in Eq. (408) by v and letting λ → λc, we have the bifurcation equation
P [u0(λc); λc]u1δu = 0, (409)
where u1 is the buckling mode deﬁned in Eq. (406).
For convenience, we introduce the following notations:
uc ≡ u0(λc) and Pc ≡ P [u0(λc); λc], (410)
then the bifurcation equation (409) adopts a more compact form
Pc u1δu = 0. (411)
8.2 Initial post-buckling analysis
Recall the principal branch u0(λ) and the bifurcated branch u0(λ) + v(λ). Let us decompose v(λ) in the
following form [67]:
v(λ) = ζu1 + ˜v(λ), (412)
94
Yang & Sharma: A tutorial on the stability and bifurcation Page 95
where u1 is the buckling mode in Eq. (406) and the bilinear inner product is u1, ˜v = 0. Then the scalar
parameter ζ can be represented by ζ = u1, v , which is a measure of the buckling mode u1 contained in
the diﬀerence v(λ) = u(λ) − u0(λ).
Similar to Eqs. (402) and (410), we introduce the following notations for convenience:
P
(n)
0 ≡ P(n)
[u0(λ); λ], P(n)
c ≡ P
(n)
0 λ=λc
,
˙P
(n)
0 ≡
∂
∂λ
P
(n)
0 ≡
∂
∂λ
P(n)
[u0(λ); λ], ˙P(n)
c ≡ ˙P
(n)
0 λ=λc
,
¨P
(n)
0 ≡
∂2
∂λ2
P
(n)
0 ≡
∂2
∂λ2
P(n)
[u0(λ); λ], ¨P(n)
c ≡ ¨P
(n)
0 λ=λc
.



(413)
Consider the Taylor series expansion of Eq. (408) at λ = λc. By Eq. (403), the expansion reads
Pc + (λ − λc) ˙Pc +
1
2
(λ − λc)2 ¨Pc + · · · vδu
+
1
2
Pc + (λ − λc) ˙Pc +
1
2
(λ − λc)2 ¨Pc + · · · v2
δu + · · · = 0.
(414)
The asymptotic expansions of the load parameter λ and ˜v(λ) in Eq. (412) for small ζ are
λ = λc + λ1ζ + λ2ζ2
+ · · · and ˜v(λ) = ζ2
u2 + ζ3
u3 + · · · , (415)
and the generalized state variable u = u0 + v = u0 + ζu1 + ˜v in Eq. (404) becomes
u = u0 + ζu1 + ζ2
u2 + ζ3
u3 + · · · . (416)
Substituting the expansions Eq. (415) into Eq. (414), we have
0 = ζ2
Pc u2 + λ1
˙Pc u1 +
1
2
Pc u2
1 δu
+ ζ3



Pc u3 + λ1
˙Pc u2 + λ2
˙Pc u1 +
1
2
λ2
1
¨Pc u1
+ Pc u1u2 +
1
2
λ1
˙Pc u2
1 +
1
6
P(4)
c u3
1



δu + o(ζ3
).
(417)
It follows that the coeﬃcients of ζ2
, ζ3
,. . . in Eq. (417) must vanish, by letting δu = u1, together
with the equality (411), we have the expressions of λ1, λ2, . . . , such that16
λ1 = −
1
2
Pc u3
1
˙Pc u2
1, (418a)
λ2 = −
1
6
P(4)
c u4
1 + Pc u2
1u2 + λ1[ ˙Pc u2u1 +
1
2
˙Pc u3
1] +
1
2
λ2
1
¨Pc u2
1
˙Pc u2
1. (418b)
In particular, for the symmetric bifurcation in many (possibly most) buckling problems [67,138], the
16See Equations (4.21) and (4.22) in the work [67]. The sign of ˙Pc u2
1 will be discussed later and the case ˙Pc u2
1 = 0 is
out of the scope of this tutorial.
95
Yang & Sharma: A tutorial on the stability and bifurcation Page 96
coeﬃcient λ1 = 0 and λ2 in Eq. (418b) reduces to
λ2 = −
1
6
P(4)
c u4
1 + Pc u2
1u2
˙Pc u2
1. (419)
8.3 Stability analysis
The diﬀerence between the energies of the state u + v and all the neighborhood states u + v + δu at a
given load parameter λ can be written as
∆P = P[u0 + v + δu; λ] − P[u0 + v; λ] =
1
2
P [u0 + v; λ](δu)2
+ o( δu 2
), (420)
where the ﬁrst-order term disappears due to equilibrium and P [u0 + v; λ](δu)2
is deﬁned in Eq. (402).
Note that the minimum of the second variation is positive, the state u + v of the system is stable at
a given load parameter λ. Thus the sign of the second variation is of primary interest. If the minimum
of the second variation is positive, then the second variation is positive for any perturbation in the
neighborhood of state u + v at λ. However, the extremum of the second variation may not exist due
to the unbound domain of the arbitrary perturbation δu. To ensure the existence of the minimum, we
discuss the minimum of the form17
P [u0 + v; λ](δu)2
δu 2
, (421)
which does not change the sign of the second variation and δu = δu, δu 1/2
.
We assume that the minimum of Eq. (421) is denoted by that exists for δu = w, such that
min
P [u0 + v; λ](δu)2
δu 2
=
P [u0 + v; λ]w2
w 2
= . (422)
Therefore, any variation w + τ ¯w admits the following inequality
P [u0 + v; λ](w + τ ¯w)2
w + τ ¯w 2
≥
P [u0 + v; λ]w2
w 2
= , (423)
alternatively,
P [u0 + v; λ]w2
+ 2τP [u0 + v; λ]w ¯w + τ2
P [u0 + v; λ] ¯w2
≥ w 2
+ 2τ w, ¯w + τ2
¯w 2
. (424)
With Eq. (422), the inequality (424) reduces to
2τ {P [u0 + v; λ]w ¯w − w, ¯w } + τ2
P [u0 + v; λ] ¯w2
− ¯w 2
≥ 0. (425)
17The extremization of this form is identical to the extremization of the second variation P [u0 + v; λ](δu)2 subjected to
the constraint (or the so-called normalization condition) δu 2 = 1. The equivalence of these two procedures can be found
in Chapter 1.18 in the book [97]. In his pioneering work [137], Koiter used the form Eq. (421) presented in this article.
96
Yang & Sharma: A tutorial on the stability and bifurcation Page 97
Since τ ∈ R is arbitrary, the inequality (425) holds if and only if
P [u0 + v; λ]w ¯w − w, ¯w = 0. (426)
For λ = λc and v = 0 and with Eq. (409) (or Eq. (411)), it is clear that w = u1 implies = 0 in Eq. (426).
Similar to the asymptotic expansions (415), in the neighborhood of λ = λc, we expand
w = u1 + ζw1 + ζ2
w2 + · · · and = ζ 1 + ζ2
2 + · · · (427)
with u1, wn = 0, n = 1, 2, 3, . . . . With the expansions in Eqs. (415) and (427), we can expand Eq. (426)
at λ = λc and u0(λc), namely
0 = ζ Pc w1 + λ1
˙Pc u1 + Pc u2
1 ¯w − 1 u1, ¯w
+ ζ2






Pc w2 + λ1
˙Pc w1 + λ2
˙Pc u1 +
1
2
λ2
1
¨Pc u1
+ Pc u1w1 + λ1
˙Pc u2
1 + Pc u2u1 +
1
2
P(4)
c u3
1


 ¯w − 1 w1, ¯w − 2 u1, ¯w



+ o(ζ2
).
(428)
Vanishing of the coeﬃcients of ζ, ζ2
, . . . gives
Pc w1 + λ1
˙Pc u1 + Pc u2
1 ¯w − 1 u1, ¯w = 0, (429a)



Pc w2 + λ1
˙Pc w1 + λ2
˙Pc u1 +
1
2
λ2
1
¨Pc u1
+ Pc u1w1 + λ1
˙Pc u2
1 + Pc u2u1 +
1
2
P(4)
c u3
1


 ¯w − 1 w1, ¯w − 2 u1, ¯w = 0. (429b)
If we choose ¯w = u1 in Eq. (429a), then we have Pc w1u1 = Pc u1w1 = 0 from Eq. (411), together
with u1 = u1, u1
1/2
= 1, we obtain
1 = λ1
˙Pc u2
1 + Pc u3
1, (430)
which, with Eq. (418a), can be further reduced to
1 = −λ1
˙Pc u2
1. (431)
Similarly, if we choose ¯w = u1 in Eq. (429b), and with Pc w2u1 = Pc u1w2 = 0 from Eq. (411) and
w1, u1 = 0 in Eq. (427), we can get
2 = λ1
˙Pc w1 + λ2
˙Pc u1 +
1
2
λ2
1
¨Pc u1 + Pc u1w1 + λ1
˙Pc u2
1 + Pc u2u1 +
1
2
P(4)
c u3
1 u1. (432)
For the case of the asymmetric bifurcation λ1 = 0, with the expansions (415) and (427) as well as the
97
Yang & Sharma: A tutorial on the stability and bifurcation Page 98
value 1 in Eq. (431), the minimum in Eq. (422) related to the second variation can be written as
= −(λ − λc) ˙Pc u2
1 + o(|λ − λc|) for asymmetric bifurcation λ1 = 0. (433)
If > 0, it corresponds to the stability of the bifurcated branch; on the contrary, the < 0 indicates
the instability of the bifurcated branch.
Again, for the case of the symmetric bifurcation λ1 = 0, it implies 1 = 0 in Eq. (431) and 2
in Eq. (432) becomes 2 = λ2
˙Pc u1 + Pc u1w1 + Pc u2u1 + 1
2 P
(4)
c u3
1 u1. Also, the zero coeﬃcient of
ζ2
in Eq. (417) gives Pc u2 + 1
2 Pc u2
1 δu = 0 while Eq. (429a) becomes Pc w1 + Pc u2
1 ¯w = 0 for
λ1 = 0 and 1 = 0. By comparing these two reduced equations and letting δu = ¯w, we then have
2u2 = w1. Now 2 further becomes 2 = λ2
˙Pc u1 + 3Pc u2u1 + 1
2 P
(4)
c u3
1 u1. Recalling the form of
λ2 in Eq. (419) for the symmetric bifurcation λ1 = 0, we ﬁnally have
2 = −2λ2
˙Pc u2
1. (434)
With Eqs. (415), (427) and (434), the minimum in Eq. (422) becomes
= −2(λ − λc) ˙Pc u2
1 + o(|λ − λc|) for symmetric bifurcation λ1 = 0. (435)
Finally, we would like to discuss the sign of ˙Pc u2
1 in Eqs. (418a), (419), (431), (433), and (435). Since
P0 u2
1 is positive for λ < λc based on the stability of the principal branch u0 below the threshold λc, and
P0 u2
1 = 0 at λ = λc from the equality (411), it follows that ˙Pc u2
1 ≤ 0. The discussion of the special case
˙Pc u2
1 = 0 is out of the scope of this tutorial and therefore the case ˙Pc u2
1 = 0 is not included here. More
discussion of the sign of ˙Pc u2
1 can be found in the work [67].
Based on the above discussion, the negative sign of ˙Pc u2
1 makes the stability conditions (433) and
(435) in alternative forms that only depend on the sign of λ1 and λ2. Dividing Eqs. (433) and (435) by
the positive value, − ˙Pc u2
1 > 0, and with the expansion of λ in Eq. (415), we have
⇒
− ˙Pc u2
1
= (λ − λc) + o(|λ − λc|) = λ1ζ + o(ζ) for asymmetric bifurcation λ1 = 0 (436)
and
⇒
− ˙Pc u2
1
= 2(λ − λc) + o(|λ − λc|) = 2λ2ζ2
+ o(ζ2
) for symmetric bifurcation λ1 = 0. (437)
Again, the value of λ1 in Eq. (436) is given by Eq. (418a) while the value of λ2 in Eq. (437) is given
by Eq. (419). Since most of the bifurcations are symmetric, then the sign of the value λ2 in Eq. (437) is
of interest. For a positive λ2, the bifurcated branch is stable. On the contrary, a negative λ2 corresponds
to an unstable bifurcated branch. The stable and unstable branches corresponding to the criteria (436)
98
Yang & Sharma: A tutorial on the stability and bifurcation Page 99
and (437) are shown in Fig. 30 in the case of a single buckling mode.
(a) (b) (c)
Figure 30: Equilibrium paths and their stabilities. Stable paths are denoted by solid curves while unstable
paths are drawn by dashed curves. The relation between the sign of the second variation and the values
λ1 and λ2 are given in Eqs. (436) and (437), respectively. (a) Asymmetric bifurcation with λ1 = 0. (b)
Symmetric bifurcation with λ2 > 0 and λ1 = 0. (c) Symmetric bifurcation with λ2 < 0 and λ1 = 0.
8.4 Example: Buckling of Euler’s column subjected to an electric ﬁeld
Let us consider the buckling and the post-buckling of an incompressible dielectric column subject to
both mechanical and the electric loads (see Fig. 31). The analysis of the purely mechanical buckling of
Euler’s column can be found in many works, to name just a few [54–56,61,65,67].
Figure 31: Schematic of the buckling of an incompressible dielectric column with pinned-pinned ends
subjected to an external electric ﬁeld E0 (in the vertical direction) and an applied force Fm
(in the
horizontal direction). The column length is L while the cross-sectional radius is r, r L. The arc
length coordinate is denoted by s. The scalar function θ(s) is the rotation (angle) of the column crosssection
from its initial horizontal direction.
99
Yang & Sharma: A tutorial on the stability and bifurcation Page 100
The potential energy of the electromechanical system shown in Fig. 31 is given by18
P[θ(s); E0, Fm
] =
L
0
1
2
B
∂θ
∂s
2
ds−
L
0
πr2
2
εE2
0 sin2
θ +
ε2
0
ε2
cos2
θ ds−Fm
L
0
(1−cos θ)ds, (438)
where θ is the rotation of the cross-section, s is the arc length coordinate, B is the bending stiﬀness, ε0
is the vacuum permittivity, ε is the material permittivity, E0 is the external electric ﬁeld, and Fm
is the
horizontal load. The state variable in Eq. (438) is regarded as the rotation θ, which can describe the
response of the column subjected to both the electric and the mechanical loads.
Consider the variation θ(s) → θ(s) + τ ¯θ(s). Here τ is a suﬃciently small parameter and ¯θ(s) is an
arbitrary function that satisﬁes all the kinematically admissible deformations. Then the ﬁrst variation
of the functional (438) is given by
δP =
d
dτ
P[θ(s) + τ ¯θ(s); E0, Fm
]
τ=0
=
L
0
B
∂θ
∂s
∂¯θ
∂s
ds −
L
0
πr2
εE2
0 1 −
ε2
0
ε2
¯θ sin θ cos θds − Fm
L
0
¯θ sin θds
= B
∂θ
∂s
¯θ
L
0
−
L
0
B
∂2
θ
∂s2
+ Fe
sin θ cos θ + Fm
sin θ ¯θds,
(439)
where Fe
is related to the electric load and deﬁned by
Fe
= πr2
εE2
0 1 −
ε2
0
ε2
. (440)
The ﬁrst variation condition, δP = 0 in Eq. (439), gives the Euler-Lagrange equation
B
∂2
θ
∂s2
+ Fe
sin θ cos θ + Fm
sin θ = 0 (441)
and the natural boundary condition
∂θ
∂s
= 0 at s = 0 & L. (442)
A trivial solution (unbuckling) to the boundary-value problem is
θ(s) = θ0(s) = 0 for 0 ≤ s ≤ L. (443)
18For the purely mechanical part, the form of the potential energy can refer to many works, e.g., [61, 65, 67]. The
second term on the RHS of Eq. (438) comes from the electric energy, which has the form − L
0 πr2 E0
0 D · dE ds for
a dielectric column subjected to an external electric ﬁeld shown in Fig. 31. The electric ﬁeld E within the column can
be decomposed into the one in the arc length and the one perpendicular to the arc length, such that E = (Es, E⊥
s ) =
(E sin θ, ε0
ε
E cos θ). Such kinds of components come from the jump conditions of the electric ﬁeld on the surface between
the vacuum and the dielectric column. For the detailed jump conditions, one can refer to the review [84]. Also, the electric
displacement D in the column has the component form D = (Ds, D⊥
s ) = (εEs, εE⊥
s ) = (εE sin θ, ε0E cos θ). It follows that
− L
0 πr2 E0
0 D · dE ds = − L
0 πr2 E0
0 εE sin2 θ +
ε2
0
ε
E cos2 θ dE ds = − L
0
πr2
2
εE2
0 sin2 θ +
ε2
0
ε2 cos2 θ ds.
100
Yang & Sharma: A tutorial on the stability and bifurcation Page 101
8.4.1 Linear bifurcation analysis
Let θ (s) be the increment of the rotation θ(s). The linearized versions of Eqs. (441) and (442) at the
trivial solution (443) are
B
∂2
θ
∂s2
+ (Fe
+ Fm
)θ = 0 (444)
and
∂θ
∂s
= 0 at s = 0 & L. (445)
Nonzero solution θ (s) of Eqs. (444) and (445) gives the condition19
Fe
+ Fm
= B
n2
π2
L2
(446)
and
θn(s) = cos
nπs
L
for 0 ≤ s ≤ L, (447)
where n ∈ Z∗
is a non-zero integer and the normalization is θn, θn ≡ 2
L
L
0
θ2
n(s)ds = 1. When
(Fe
+ Fm
) increases from zero to Bπ2
/L2
, the ﬁrst buckling mode θ1(s) = cos πs
L occurs.
8.4.2 Energy stability criterion
In contrast to the linear bifurcation analysis, we study the stability of the electromechanical system
by the energy method. Thus we have to examine the second variation condition. Similar to the ﬁrst
variation (439), the second variation of the potential energy (438) is
δ2
P =
d2
dτ2
P[θ(s) + τ ¯θ(s); E0, Fm
]
τ=0
=
L
0
B
∂¯θ
∂s
2
ds − Fe
L
0
¯θ2
cos(2θ)ds − Fm
L
0
¯θ2
cos θds.
(448)
At the trivial solution θ0(s) = 0 in Eq. (443), the second variation (448) reads
δ2
P|θ0=0 =
L
0
B
∂¯θ
∂s
2
− (Fe
+ Fm
)¯θ2
ds. (449)
The integral (449) is a continuous functional of the smooth variation ¯θ(s). If the minimum value of
the functional (449) is non-negative, then the equilibrium solution θ0(s) = 0 in Eq. (443) is stable. To
ensure the existence of the minimum value of the second variation δ2
P|θ0=0 in Eq. (449), we introduce
a compact subset of the domain of the functional, which does not change the positive deﬁniteness of the
functional. Consider the following normalization condition
¯θ 2
= ¯θ, ¯θ ≡
2
L
L
0
¯θ2
ds = 1. (450)
Based on the boundary condition (442), we consider the Fourier cosine series of the kinematically ad-
19Note that the linear bifurcation analysis only gives the necessary condition for the existence of non-trivial solutions.
101
Yang & Sharma: A tutorial on the stability and bifurcation Page 102
missible variation
¯θ(s) =
∞
n=0
an cos
nπs
L
. (451)
The normalization condition (450), with the orthogonality condition, now reads
∞
n=0
a2
n = 1. (452)
With Eqs. (450)−(452), the integral (449) is bounded and becomes
δ2
P|θ0=0 =
L
2
∞
n=0
a2
n B
nπ
L
2
− (Fe
+ Fm
) ≥ , (453)
where is the minimum of the integral, namely
= min δ2
P|θ0=0 =
L
2
B
π
L
2
− (Fe
+ Fm
) at ¯θ = cos
πs
L
. (454)
Apparently, the minimum value in Eq. (454) is positive for 0 < (Fe
+ Fm
) < B π
L
2
, which
indicates that the trivial solution θ0(s) = 0 is stable. On the contrary, the trivial solution θ0(s) = 0
becomes unstable for (Fe
+ Fm
) > B π
L
2
since the minimum is negative.
8.4.3 Post-buckling analysis: LSK approach
By the notations in Eq. (402), the ﬁrst variation of the potential energy (438) is
P [θ; E0, Fm
]θ1 =
L
0
B
∂θ
∂s
∂θ1
∂s
ds − Fe
L
0
θ1 sin θ cos θds − Fm
L
0
θ1 sin θds, (455)
and the higher-order variations are
P [θ; E0, Fm
]θ1θ2 =
L
0
B
∂θ2
∂s
∂θ1
∂s
ds − Fe
L
0
θ1θ2 cos(2θ)ds − Fm
L
0
θ1θ2 cos θds, (456)
P [θ; E0, Fm
]θ1θ2θ3 = 2Fe
L
0
θ1θ2θ3 sin(2θ)ds + Fm
L
0
θ1θ2θ3 sin θds, (457)
P(4)
[θ; E0, Fm
]θ1θ2θ3θ4 = 4Fe
L
0
θ1θ2θ3θ4 cos(2θ)ds + Fm
L
0
θ1θ2θ3θ4 cos θds. (458)
Note that the generalized displacement u in Eq. (402) is the rotation θ(s) while the load parameter λ
is the horizontal load Fm
. The electric ﬁeld E0 is regarded as a constant in this example for simplicity.
Thus we only have one control parameter Fm
here.
Consider the trivial solution θ(s) = θ0(s) = 0, i.e., the prebuckling solution of the column. With the
second variation (456), the bifurcation equation (409) becomes
0 = P [0; E0, Fm
c ]θ1δθ =
L
0
B
∂δθ
∂s
∂θ1
∂s
ds − Fe
L
0
θ1δθds − Fm
c
L
0
θ1δθds, (459)
102
Yang & Sharma: A tutorial on the stability and bifurcation Page 103
which gives similar results as that of the linear bifurcation analysis, namely
B
∂2
θ1
∂s2
+ (Fe
+ Fm
c )θ1 = 0 (460)
and the natural boundary condition
∂θ1
∂s
= 0 at s = 0 & L. (461)
The lowest eigenvalue corresponds to the critical load, namely20
Fe
+ Fm
c =
Bπ2
L2
, (462)
and the normalized eigenfunction θ1(s), with the normalization θ1, θ1 ≡ 2
L
L
0
θ2
1ds = 1 in Eq. (450), is
θ1 = cos
πs
L
. (463)
With the third and fourth variations (457) and (458) as well as Eq. (463), we obtain
Pc θ3
1 = P [0; E0, Fm
c ]θ3
1 = 0, (464a)
P(4)
c θ4
1 = P(4)
[0; E0, Fm
c ]θ4
1 = (4Fe
+ Fm
c )
L
0
θ4
1ds =
3
8
(4Fe
+ Fm
c )L. (464b)
By Eq. (457) and θ(s) = θ0(s) = 0 for the principal branch (unbuckled), we have
Pc θ2
1θ2 = P [0; E0, Fm
c ]θ2
1θ2 = 0. (464c)
Also, with Eqs. (456) and (463), we have the partial derivative
˙Pc θ2
1 =
∂
∂Fm
P [θ]θ2
1
(θ0,F m
c )
= −
L
0
θ2
1ds = −
L
2
. (464d)
Now from Eq. (464), we can get the values of λ1 in Eq. (418a) and λ2 in Eq. (419), such that
λ1 = −
1
2
Pc θ3
1
˙Pc θ2
1 = 0, (465a)
λ2 = −
1
6
P(4)
c θ4
1 + Pc θ2
1θ2
˙Pc θ2
1 =
1
8
(4Fe
+ Fm
c ). (465b)
Since λ1 = 0 and λ2 > 0, from Eq. (437) and the schematic in Fig. 30(b), the bifurcation branch is
symmetric and stable after the horizontal load Fm
slightly exceeds the threshold Fm
c . With Eq. (465),
20For simplicity, we just consider the case in which the electric load Fe is 0 < Fe < Bπ2
L2 . And there is no buckling if
the applied mechanical load Fm is below the critical load Fm
c .
103
Yang & Sharma: A tutorial on the stability and bifurcation Page 104
the expansion of the load parameter Fm
is
Fm
= Fm
c + λ1ζ + λ2ζ2
+ o(ζ2
) = Fm
c +
1
8
(4Fe
+ Fm
c )ζ2
+ o(ζ2
).
Since the small scalar ζ is the amplitude of the buckling mode θ1 in Eq. (463), we have
Fm
Fm
c
= 1 +
1
8
1 + 4
Fe
Fm
c
ζ2
(466)
and the load-shortening distance is
L =
L
0
1 − cos ζθ1(s) ds ∼=
L
0
ζ2
θ2
1(s)
2
ds =
ζ2
L
4
. (467)
By Eq. (467), Eq. (466) can be rewritten as
Fm
Fm
c
= 1 +
1
2
L
L
1 + 4
Fe
Fm
c
. (468)
From Fe
+Fm
c = Bπ2
L2 in Eq. (462), we introduce the ratio β = Fe Bπ2
L2 , and then Fm
c = (1−β)Bπ2
L2 .
It follows that
Fm
Fm
c
= 1 +
1
2
L
L
1 + 4
β
1 − β
. (469)
For β = 0, we recover the classic result for the post-buckling analysis of Euler’s column subjected to
purely mechanical loading [61,65].
9 Concluding remarks and future study
In what is intended to be a simple tutorial article, we can hope to be neither comprehensive nor as
rigorous as we would like. As highlighted earlier in the introduction, the reader should take the current
tutorial merely as a starting point before delving into the rather rich literature that already exists on
the topic of instability, bifurcation, and electromechanical coupling. Further study of these topics could
be taken in several directions which we outline below along with some selected references.
While the ﬁrst section of our tutorial presented a broad overview of the ﬁeld of stability and bifurcation
analysis, our primary focus, and all the relevant examples, are in the context of elastic materials
that do not exhibit dissipation. In soft materials that exhibit some sort of electromechanical coupling,
dissipation can emanate from diﬀusion of ions, electronic conduction, formation of defects, or viscous
eﬀects in the mechanical behavior of the polymer. Invoking viscoelasticity, current leakage and dielectric
relaxation, stability analysis of dielectric elastomers was conducted by Chiang et. al. [147, 148]. A detailed
examination of the eﬀect of viscoelasticity on electromechanical instabilities of dielectric elastomers
was conducted by Wang et al. [149] where a ﬁnite element implementation was also developed.
A comprehensive theoretical framework for dissipation due to diﬀusing ions in dielectrics was created
104
Yang & Sharma: A tutorial on the stability and bifurcation Page 105
by Xiao and Bhattacharya [150] while Darbaniyan et. al. [151] provide a model with thermal eﬀects.
Although other works exist this topic, these two aforementioned references (and citations therein) may
be used to acquire a broad overview of modeling eﬀorts in diﬀusion and thermal related dissipation.
However, neither of these two works focus on instabilities and examination of thermal and diﬀusive instabilities
in electromechanical soft materials remains a (relatively) understudied sub-ﬁeld.
Another important aspect pertaining to instabilities in soft matter is related to extreme deformations
where a singularity may exist, for example, electro-creasing instability involves a singularity. Similar
to [41], electro-creasing could have applications in capacitors and energy harvesting; however, an analytical
solution of creasing instability is still an open problem. A related example is the problem related
to the crumpling of thin sheets [152–155] that may be harnessed for energy harvesting in wearable elec-
tronics.
A signiﬁcant line of research is ongoing on developing soft composites that exhibit unusual macroscopic
electromechanical behavior. Speciﬁcally, germane to the current article, several researchers have
investigated how instabilities at the microscale may inﬂuence the overall response of the heterogeneous
materials. For example, Ponte Casta˜neda and Siboni [156] proposed a homogenization framework for
electro-elastic composite materials at ﬁnite strains incorporating nonlinear dielectric behavior. Using
this homogenization framework, electromechanical instabilities in ﬁber-constrained, dielectric-elastomer
composites subjected to electromechanical loadings were further studied [157–159]. In addition, electromechanical
instabilities in soft heterogeneous dielectric elastomers were investigated in multi-layered
dielectric composites at ﬁnite strains, eg., Bertoldi and Gei [103], Rudykh and deBotton [160], Rudykh
et al. [161], Spinelli and Lopez-Pamies [162].
We remark that the solution for most nonlinear problems in electroelasticity, and certainly that involve
stability and bifurcation analysis, necessarily require computational approaches. In our tutorial, in
order to highlight the basic concepts, we have restricted ourselves to examples that are amenable to analytical
solutions. In order to explore the computational literature, we refer the reader to Vu et al. [163]
who presented a ﬁnite element approach to simulate the nonlinear coupling behavior of electroactive
polymers under electric stimulation. To investigate the inﬂuence of the free space on the electric ﬁeld
and on the deformation ﬁeld, Vu and Steinmann [164] employed a coupled BEM-FEM approach to exploit
the advantage of the ﬁnite element method in solving nonlinear problems in nonlinear electroelasticity.
For dissipative electro-magneto-mechanics, Miehe et al. [165], based on incremental variational principles,
presented a general framework for the macroscopic, continuum-based formulation and numerical
implementation of dissipative functional materials with electro-magneto-mechanical couplings. For dissipative
electro-mechanics, Z¨ah and Miehe [166] developed an approach for computational homogenization
based on rigorous exploitation of rate-type and incremental variational principles. In ﬁnite magnetoelectro-elasticity,
Miehe et al. [167,168] outlined a variational-based framework for homogenization and
multiscale stability analysis, which allows tracking of post-critical solution paths such as those related to
105
Yang & Sharma: A tutorial on the stability and bifurcation Page 106
pull-in instabilities. Polukhov et al. [169] investigated the computational stability analysis of electroactive
polymer composites by implementing Bloch-Floquet wave analysis in the context of a ﬁnite element
discretization, with a particular focus on the investigation of macroscopic loss of strong ellipticity and
microscopic bifurcation-type instabilities.
The conﬂuence of stability and electromechanical eﬀects plays a signiﬁcant role in biological cells
also—arguably an important subset of soft matter. Among other consequences of this, such an interaction
between cells and electric ﬁelds involves changes in the shapes of the cells as well as formation
of pores (i.e., electroporation). Aside from the classical works by Helfrich [170, 171], the following are
some recent representative works on this topic which also contain citations to the vast literature on this
subject [172–174].
Finally, we remark that despite numerous works in the purely mechanical case, very few have addressed
post-buckling of soft dielectrics.
Acknowledgments
Support from the University of Houston and the M. D. Anderson Professorship is gratefully acknowledged.
S.Y. would like to acknowledge Qilu Youth Scholar Program of Shandong University and the Natural
Science Foundation of Shandong Province of China (Grant No. ZR2017MA037).
References
[1] F. Ilievski, A. D. Mazzeo, R. F. Shepherd, X. Chen, and G. M. Whitesides, “Soft robotics for
chemists,” Angewandte Chemie, vol. 123, no. 8, pp. 1930–1935, 2011.
[2] D. Yang, B. Mosadegh, A. Ainla, B. Lee, F. Khashai, Z. Suo, K. Bertoldi, and G. M. Whitesides,
“Buckling of elastomeric beams enables actuation of soft machines,” Advanced Materials, vol. 27,
no. 41, pp. 6323–6327, 2015.
[3] S. Bauer, S. Bauer-Gogonea, I. Graz, M. Kaltenbrunner, C. Keplinger, and R. Schw¨odiauer, “25th
anniversary article: a soft future: from robots and sensor skin to energy harvesters,” Advanced
Materials, vol. 26, no. 1, pp. 149–162, 2014.
[4] R. Pelrine, R. D. Kornbluh, J. Eckerle, P. Jeuck, S. Oh, Q. Pei, and S. Stanford, “Dielectric
elastomers: generator mode fundamentals and applications,” in Smart Structures and Materials
2001: Electroactive Polymer Actuators and Devices (Y. Bar-Cohen, ed.), vol. 4329, pp. 148 – 156,
International Society for Optics and Photonics, SPIE, 2001.
[5] S. J. A. Koh, X. Zhao, and Z. Suo, “Maximal energy that can be converted by a dielectric elastomer
generator,” Applied Physics Letters, vol. 94, no. 26, p. 262902, 2009.
106
Yang & Sharma: A tutorial on the stability and bifurcation Page 107
[6] S. J. A. Koh, C. Keplinger, T. Li, S. Bauer, and Z. Suo, “Dielectric elastomer generators: How much
energy can be converted?,” IEEE/ASME Transactions on mechatronics, vol. 16, no. 1, pp. 33–41,
2011.
[7] F. Invernizzi, S. Dulio, M. Patrini, G. Guizzetti, and P. Mustarelli, “Energy harvesting from human
motion: materials and techniques,” Chem. Soc. Rev., vol. 45, pp. 5455–5473, 2016.
[8] G. Kofod, W. Wirges, M. Paajanen, and S. Bauer, “Energy minimization for self-organized structure
formation and actuation,” Applied Physics Letters, vol. 90, no. 8, p. 081916, 2007.
[9] F. Carpi, S. Bauer, and D. De Rossi, “Stretching dielectric elastomer performance,” Science,
vol. 330, no. 6012, pp. 1759–1761, 2010.
[10] C. Keplinger, M. Kaltenbrunner, N. Arnold, and S. Bauer, “R¨ontgens electrode-free elastomer
actuators without electromechanical pull-in instability,” Proceedings of the National Academy of
Sciences, vol. 107, no. 10, pp. 4505–4510, 2010.
[11] H. Shao, S. Wei, X. Jiang, D. P. Holmes, and T. K. Ghosh, “Bioinspired electrically activated soft
bistable actuators,” Advanced Functional Materials, vol. 28, no. 35, p. 1802999, 2018.
[12] D. van den Ende, J.-D. Kamminga, A. Boersma, T. Andritsch, and P. G. Steeneken, “Voltagecontrolled
surface wrinkling of elastomeric coatings,” Advanced Materials, vol. 25, no. 25, pp. 3438–
3442, 2013.
[13] R. Pelrine, R. Kornbluh, Q. Pei, and J. Joseph, “High-speed electrically actuated elastomers with
strain greater than 100%,” Science, vol. 287, no. 5454, pp. 836–839, 2000.
[14] J. Huang, T. Li, C. Chiang Foo, J. Zhu, D. R. Clarke, and Z. Suo, “Giant, voltage-actuated
deformation of a dielectric elastomer under dead load,” Applied Physics Letters, vol. 100, no. 4,
p. 041911, 2012.
[15] X. Zhao and Z. Suo, “Theory of dielectric elastomers capable of giant deformation of actuation,”
Physical review letters, vol. 104, no. 17, p. 178302, 2010.
[16] M. A. Biot, “Surface instability of rubber in compression,” Applied Scientiﬁc Research, Section A,
vol. 12, no. 2, pp. 168–182, 1963.
[17] N. Bowden, S. Brittain, A. G. Evans, J. W. Hutchinson, and G. M. Whitesides, “Spontaneous
formation of ordered structures in thin ﬁlms of metals supported on an elastomeric polymer,”
Nature, vol. 393, no. 6681, pp. 146–149, 1998.
[18] S. Yang, K. Khare, and P.-C. Lin, “Harnessing surface wrinkle patterns in soft matter,” Advanced
Functional Materials, vol. 20, no. 16, pp. 2550–2564, 2010.
[19] Q. Wang and X. Zhao, “Creasing-wrinkling transition in elastomer ﬁlms under electric ﬁelds,”
Physical Review E, vol. 88, no. 4, p. 042403, 2013.
107
Yang & Sharma: A tutorial on the stability and bifurcation Page 108
[20] X. Liu, B. Li, H. Chen, S. Jia, and J. Zhou, “Voltage-induced wrinkling behavior of dielectric
elastomer,” Journal of Applied Polymer Science, vol. 133, no. 14, 2016.
[21] H. Godaba, Z.-Q. Zhang, U. Gupta, C. C. Foo, and J. Zhu, “Dynamic pattern of wrinkles in a
dielectric elastomer,” Soft matter, vol. 13, no. 16, pp. 2942–2951, 2017.
[22] A. Dorfmann and R. W. Ogden, “Nonlinear electroelastostatics: Incremental equations and stability,”
International Journal of Engineering Science, vol. 48, no. 1, pp. 1–14, 2010.
[23] A. Gent and I. Cho, “Surface instabilities in compressed or bent rubber blocks,” Rubber Chemistry
and Technology, vol. 72, no. 2, pp. 253–262, 1999.
[24] V. Trujillo, J. Kim, and R. C. Hayward, “Creasing instability of surface-attached hydrogels,” Soft
Matter, vol. 4, no. 3, pp. 564–569, 2008.
[25] W. Hong, X. Zhao, and Z. Suo, “Formation of creases on the surfaces of elastomers and gels,”
Applied Physics Letters, vol. 95, no. 11, p. 111901, 2009.
[26] E. Hohlfeld and L. Mahadevan, “Unfolding the sulcus,” Physical review letters, vol. 106, no. 10,
p. 105702, 2011.
[27] Q. Wang, L. Zhang, and X. Zhao, “Creasing to cratering instability in polymers under ultrahigh
electric ﬁelds,” Physical review letters, vol. 106, no. 11, p. 118301, 2011.
[28] Q. Wang, M. Tahir, J. Zang, and X. Zhao, “Dynamic electrostatic lithography: Multiscale ondemand
patterning on large-area curved surfaces,” Advanced Materials, vol. 24, no. 15, pp. 1947–
1951, 2012.
[29] G. Zurlo, M. Destrade, D. DeTommasi, and G. Puglisi, “Catastrophic thinning of dielectric elastomers,”
Physical Review Letters, vol. 118, no. 7, p. 078001, 2017.
[30] L. Pocivavsek, R. Dellsy, A. Kern, S. Johnson, B. Lin, K. Y. C. Lee, and E. Cerda, “Stress and
fold localization in thin elastic membranes,” Science, vol. 320, no. 5878, pp. 912–916, 2008.
[31] P. Kim, M. Abkarian, and H. A. Stone, “Hierarchical folding of elastic membranes under biaxial
compressive stress,” Nature materials, vol. 10, no. 12, pp. 952–957, 2011.
[32] B. Tavakol, M. Bozlar, C. Punckt, G. Froehlicher, H. A. Stone, I. A. Aksay, and D. P. Holmes,
“Buckling of dielectric elastomeric plates for soft, electrically active microﬂuidic pumps,” Soft
Matter, vol. 10, no. 27, pp. 4789–4794, 2014.
[33] B. Tavakol and D. P. Holmes, “Voltage-induced buckling of dielectric ﬁlms using ﬂuid electrodes,”
Applied Physics Letters, vol. 108, no. 11, p. 112901, 2016.
[34] H. Bense, M. Trejo, E. Reyssat, J. Bico, and B. Roman, “Buckling of elastomer sheets under
non-uniform electro-actuation,” Soft matter, vol. 13, no. 15, pp. 2876–2885, 2017.
108
Yang & Sharma: A tutorial on the stability and bifurcation Page 109
[35] S. Yang, X. Zhao, and P. Sharma, “Revisiting the instability and bifurcation behavior of soft
dielectrics,” Journal of Applied Mechanics-Transactions of the ASME, vol. 84, p. 031008, MARCH
2017.
[36] K. Stark and C. Garton, “Electric strength of irradiated polythene,” Nature, vol. 176, no. 4495,
pp. 1225–1226, 1955.
[37] J.-S. Plante and S. Dubowsky, “Large-scale failure modes of dielectric elastomer actuators,” International
journal of solids and structures, vol. 43, no. 25, pp. 7727–7751, 2006.
[38] X. Zhao and Z. Suo, “Method to analyze electromechanical stability of dielectric elastomers,”
Applied Physics Letters, vol. 91, no. 6, p. 061921, 2007.
[39] X. Zhao and Z. Suo, “Electromechanical instability in semicrystalline polymers,” Applied Physics
Letters, vol. 95, no. 3, p. 031904, 2009.
[40] S. Yang, X. Zhao, and P. Sharma, “Avoiding the pull-in instability of a dielectric elastomer ﬁlm and
the potential for increased actuation and energy harvesting,” Soft Matter, vol. 13, pp. 4552–4558,
2017.
[41] Q. Wang, Z. Suo, and X. Zhao, “Bursting drops in solid dielectrics caused by high voltages,” Nature
communications, vol. 3, p. 1157, 2012.
[42] X. Zhao and Q. Wang, “Harnessing large deformation and instabilities of soft dielectrics: theory,
experiment, and application,” Applied Physics Reviews, vol. 1, no. 2, p. 021304, 2014.
[43] N. Hu and R. Burgue˜no, “Buckling-induced smart applications: recent advances and trends,”
Smart Materials and Structures, vol. 24, no. 6, p. 063001, 2015.
[44] K. Bertoldi, P. M. Reis, S. Willshaw, and T. Mullin, “Negative poisson’s ratio behavior induced
by an elastic instability,” Advanced Materials, vol. 22, no. 3, pp. 361–366, 2010.
[45] E. P. Chan, E. J. Smith, R. C. Hayward, and A. J. Crosby, “Surface wrinkles for smart adhesion,”
Advanced Materials, vol. 20, no. 4, pp. 711–716, 2008.
[46] Y. Forterre, J. M. Skotheim, J. Dumais, and L. Mahadevan, “How the venus ﬂytrap snaps,” Nature,
vol. 433, no. 7024, p. 421, 2005.
[47] D. P. Holmes, “Elasticity and stability of shape-shifting structures,” Current Opinion in Colloid
& Interface Science, vol. 40, pp. 118 – 137, 2019. Particle Systems.
[48] S. H. Strogatz, Nonlinear dynamics and chaos: with applications to physics, biology, chemistry,
and engineering. Westview press, 2014.
[49] R. Seydel, Practical bifurcation and stability analysis, vol. 5. Springer Science & Business Media,
2009.
109
Yang & Sharma: A tutorial on the stability and bifurcation Page 110
[50] J. Guckenheimer and P. Holmes, Nonlinear oscillations, dynamical systems, and bifurcations of
vector ﬁelds, vol. 42. Springer Science & Business Media, 2013.
[51] L. Perko, Diﬀerential equations and dynamical systems, vol. 7. Springer Science & Business Media,
2013.
[52] Y. A. Kuznetsov, Elements of applied bifurcation theory, vol. 112. Springer Science & Business
Media, 2013.
[53] S. Sastry, Nonlinear systems: analysis, stability, and control, vol. 10. Springer Science & Business
Media, 2013.
[54] S. P. Timoshenko and J. M. Gere, Theory of Elastic Stability. McGrawHill-Kogakusha Ltd, Tokyo,
1961.
[55] Z. P. BAˇZANT and L. CEDOLIN, Stability of structures: elastic, inelastic, fracture and damage
theories. World Scientiﬁc, 2010.
[56] S. Antman, Nonlinear Problems of Elasticity, vol. 107. Springer Science & Business Media, 2013.
[57] R. W. Ogden, Non-linear elastic deformations. Courier Corporation, 1997.
[58] L. Dorfmann and R. W. Ogden, Nonlinear theory of electroelastic and magnetoelastic interactions.
Springer Science & Business Media, 2014.
[59] L. Dorfmann and R. W. Ogden, “Instabilities of soft dielectrics,” Philosophical Transactions of the
Royal Society A: Mathematical, Physical and Engineering Sciences, vol. 377, no. 2144, p. 20180077,
2019.
[60] M. Golubitsky, I. Stewart, and D. G. Schaeﬀer, Singularities and groups in bifurcation theory,
vol. 2. Springer Science & Business Media, 2012.
[61] A. M. Van der Heijden, WT Koiter’s elastic stability of solids and structures, vol. 10. Cambridge
University Press Cambridge, 2009.
[62] L. Liu, “On energy formulations of electrostatics for continuum media,” Journal of the Mechanics
and Physics of Solids, vol. 61, no. 4, pp. 968–990, 2013.
[63] L. Liu, “An energy formulation of continuum magneto-electro-elasticity with applications,” Journal
of the Mechanics and Physics of Solids, vol. 63, pp. 451–480, 2014.
[64] A. M. Lyapunov, “The general problem of the stability of motion,” International journal of control,
vol. 55, no. 3, pp. 531–773, 1992.
[65] N. Triantafyllidis, Stability of solids: From structures to materials. Ecole polytechnique, 2011.
[66] J. M. T. Thompson and G. W. Hunt, A general theory of elastic stability. Wiley, 1973.
110
Yang & Sharma: A tutorial on the stability and bifurcation Page 111
[67] B. Budiansky, “Theory of buckling and post-buckling behavior of elastic structures,” in Advances
in applied mechanics, vol. 14, pp. 1–65, Elsevier, 1974.
[68] R. J. Knops and E. Wilkes, “Theory of elastic stability(liapunov functions application to stability
analysis of dynamic systems and elastic bodies, considering eigenfunction method, maximum principle
and energy criterion),” Solid-state mechanics 3.(A 73-45495 24-32) Berlin, Springer-Verlag,
1973,, pp. 125–302, 1973.
[69] M. Como and A. Grimaldi, Theory of stability of continuous elastic structures, vol. 1. CRC Press,
1995.
[70] Y.-c. Chen, S. Yang, and L. Wheeler, “Surface instability of elastic half-spaces by using the energy
method,” Proc. R. Soc. A, vol. 474, no. 2213, p. 20170854, 2018.
[71] R. Courant and D. Hilbert, Methods of mathematical physics, Volume I. Interscience Publishers,
1953.
[72] R. Weinstock, Calculus of variations: with applications to physics and engineering. Courier Corporation,
1974.
[73] J. M. Ball, “The calculus of variations and materials science,” Quarterly of Applied Mathematics,
vol. 56, no. 4, pp. 719–740, 1998.
[74] Y. Chen, Nonlinear Elasticity: Theory and Applications, ch. Singularity Theory and Nonlinear
Bifurcation Analysis, pp. 305–344. Cambridge University Press, Cambridge. UK., 2001.
[75] J. Ericksen and R. Toupin, “Implications of hadamard’s conditions for elastic stability with respect
to uniqueness theorems,” Canadian Journal of Mathematics, vol. 8, pp. 432–436, 1956.
[76] R. Hill, “On uniqueness and stability in the theory of ﬁnite elastic strain,” Journal of the Mechanics
and Physics of Solids, vol. 5, no. 4, pp. 229–241, 1957.
[77] Y.-C. Chen and D. Haughton, “Stability and bifurcation of inﬂation of elastic cylinders,” in Proceedings
of the Royal Society of London A: Mathematical, Physical and Engineering Sciences, vol. 459,
pp. 137–156, The Royal Society, 2003.
[78] Y.-c. Chen and E. Fried, “Stability and bifurcation of a soap ﬁlm spanning a ﬂexible loop,” Journal
of Elasticity, vol. 116, no. 1, pp. 75–100, 2014.
[79] M. E. Gurtin, An introduction to continuum mechanics, vol. 158. Academic press, 1982.
[80] M. E. Gurtin, E. Fried, and L. Anand, The mechanics and thermodynamics of continua. Cambridge
University Press, 2010.
[81] A. G. Holzapfel, Nonlinear Solid Mechanics II. John Wiley & Sons, Inc., 2000.
[82] R. M. McMeeking and C. M. Landis, “Electrostatic forces and stored energy for deformable dielectric
materials,” Journal of Applied Mechanics, vol. 72, no. 4, pp. 581–590, 2005.
111
Yang & Sharma: A tutorial on the stability and bifurcation Page 112
[83] R. Bustamante, A. Dorfmann, and R. W. Ogden, “On electric body forces and maxwell stresses in
nonlinearly electroelastic solids,” International Journal of Engineering Science, vol. 47, no. 11-12,
pp. 1131–1141, 2009.
[84] L. Dorfmann and R. W. Ogden, “Nonlinear electroelasticity: material properties, continuum theory
and applications,” in Proc. R. Soc. A, vol. 473, p. 20170311, The Royal Society, 2017.
[85] C. Truesdell and W. Noll, “The non-linear ﬁeld theories of mechanics,” in The non-linear ﬁeld
theories of mechanics, pp. 1–579, Springer, 2004.
[86] J. N. Reddy, An introduction to continuum mechanics. Cambridge university press, 2007.
[87] Y. C. Fung and P. Tong, Classical and computational solid mechanics, vol. 1. World Scientiﬁc
Publishing Co Inc, 2001.
[88] Z. Huang, Fundamentals of continuum mechanics. Higher Education Press, Beijing (in Chinese),
2003.
[89] R. A. Toupin, “The elastic dielectric,” Journal of Rational Mechanics and Analysis, vol. 5, no. 6,
pp. 849–915, 1956.
[90] R. Toupin, “Stress tensors in elastic dielectrics,” Archive for Rational Mechanics and Analysis,
vol. 5, no. 1, p. 440, 1960.
[91] J. Ericksen, “Electromagnetic eﬀects in thermoelastic materials,” Mathematics and mechanics of
solids, vol. 7, no. 2, pp. 165–189, 2002.
[92] A. Dorfmann and R. Ogden, “Nonlinear electroelasticity,” Acta Mechanica, vol. 174, no. 3-4,
pp. 167–183, 2005.
[93] J. Ericksen, “Theory of elastic dielectrics revisited,” Archive for Rational Mechanics and Analysis,
vol. 183, no. 2, pp. 299–313, 2007.
[94] R. Fosdick and H. Tang, “Electrodynamics and thermomechanics of material bodies,” Journal of
Elasticity, vol. 88, no. 3, pp. 255–297, 2007.
[95] Z. Suo, X. Zhao, and W. H. Greene, “A nonlinear ﬁeld theory of deformable dielectrics,” Journal
of the Mechanics and Physics of Solids, vol. 56, no. 2, pp. 467–486, 2008.
[96] Z. Suo, “Theory of dielectric elastomers,” Acta Mechanica Solida Sinica, vol. 23, no. 6, pp. 549–578,
2010.
[97] K. F. Riley and M. P. Hobson, Essential mathematical methods for the physical sciences. Cambridge
University Press, 2011.
[98] C. D. Meyer, Matrix analysis and applied linear algebra, vol. 2. Siam, 2000.
[99] M. Mooney, “A theory of large elastic deformation,” Journal of applied physics, vol. 11, no. 9,
pp. 582–592, 1940.
112
Yang & Sharma: A tutorial on the stability and bifurcation Page 113
[100] R. Rivlin, “Large elastic deformations of isotropic materials. iv. further developments of the general
theory,” Philosophical Transactions of the Royal Society of London A: Mathematical, Physical and
Engineering Sciences, vol. 241, no. 835, pp. 379–397, 1948.
[101] A. Gent, “A new constitutive relation for rubber,” Rubber chemistry and technology, vol. 69, no. 1,
pp. 59–61, 1996.
[102] X. Zhao, W. Hong, and Z. Suo, “Electromechanical hysteresis and coexistent states in dielectric
elastomers,” Physical review B, vol. 76, no. 13, p. 134113, 2007.
[103] K. Bertoldi and M. Gei, “Instabilities in multilayered soft dielectrics,” Journal of the Mechanics
and Physics of Solids, vol. 59, no. 1, pp. 18 – 42, 2011.
[104] H. S. Park, Q. Wang, X. Zhao, and P. A. Klein, “Electromechanical instability on dielectric polymer
surface: Modeling and experiment,” Computer Methods in Applied Mechanics and Engineering,
vol. 260, pp. 40–49, 2013.
[105] L. Dorfmann and R. W. Ogden, “Instabilities of an electroelastic plate,” International Journal of
Engineering Science, vol. 77, pp. 79–101, 2014.
[106] T. Lu, L. An, J. Li, C. Yuan, and T. Wang, “Electro-mechanical coupling bifurcation and bulging
propagation in a cylindrical dielectric elastomer tube,” Journal of the Mechanics and Physics of
Solids, vol. 85, pp. 160–175, 2015.
[107] X. Liang and S. Cai, “New electromechanical instability modes in dielectric elastomer balloons,”
International Journal of Solids and Structures, 2017.
[108] S. Seiﬁ and H. Park, “Electro-elastocapillary rayleigh-plateau instability in dielectric elastomer
ﬁlms,” Soft Matter, vol. 13, pp. 4305–4310, 2017.
[109] E. Bortot and G. Shmuel, “Tuning sound with soft dielectrics,” Smart Materials and Structures,
vol. 26, p. 045028, 04 2017.
[110] E. Bortot and G. Shmuel, “Prismatic bifurcations of soft dielectric tubes,” International Journal
of Engineering Science, vol. 124, pp. 104–114, 2018.
[111] G. Mao, L. Wu, X. Liang, and S. Qu, “Morphology of voltage-triggered ordered wrinkles of a
dielectric elastomer sheet,” Journal of Applied Mechanics, vol. 84, no. 11, p. 111005, 2017.
[112] G. Mao, L. Wu, Y. Fu, J. Liu, and S. Qu, “Voltage-controlled radial wrinkles of a trumpet-like
dielectric elastomer structure,” AIP Advances, vol. 8, no. 3, p. 035314, 2018.
[113] Y. Su, H. C. Broderick, W. Chen, and M. Destrade, “Wrinkles in soft dielectric plates,” Journal
of the Mechanics and Physics of Solids, vol. 119, pp. 298–318, 2018.
[114] Y. Su, B. Wu, W. Chen, and M. Destrade, “Finite bending and pattern evolution of the associated
instability for a dielectric elastomer slab,” International Journal of Solids and Structures, vol. 158,
pp. 191–209, 2019.
113
Yang & Sharma: A tutorial on the stability and bifurcation Page 114
[115] P. Greaney, M. Meere, and G. Zurlo, “The out-of-plane behaviour of dielectric membranes: description
of wrinkling and pull-in instabilities,” Journal of the Mechanics and Physics of Solids,
vol. 122, pp. 84–97, 2019.
[116] Y. Fu, Y. Xie, and L. Dorfmann, “A reduced model for electrodes-coated dielectric plates,” International
Journal of Non-Linear Mechanics, 2018.
[117] R. Huang and Z. Suo, “Electromechanical phase transition in dielectric elastomers,” vol. 468,
pp. 1014–1040, The Royal Society, 2012.
[118] T. He, X. Zhao, and Z. Suo, “Dielectric elastomer membranes undergoing inhomogeneous deformation,”
Journal of Applied Physics, vol. 106, no. 8, p. 083522, 2009.
[119] S. Yang and Y.-c. Chen, “Wrinkle surface instability of an inhomogeneous elastic block with graded
stiﬀness,” in Proc. R. Soc. A, vol. 473, p. 20160882, The Royal Society, 2017.
[120] A. Norris, “Comment on“method to analyze electromechanical stability of dielectric elastomers”[appl.
phys. lett. 91, 061921 (2007)],” Applied Physics Letters, vol. 92, no. 2, p. 026101,
2008.
[121] R. D´ıaz-Calleja, E. Riande, and M. Sanchis, “On electromechanical stability of dielectric elastomers,”
Applied Physics Letters, vol. 93, no. 10, p. 101902, 2008.
[122] B.-X. Xu, R. Mueller, M. Klassen, and D. Gross, “On electromechanical stability analysis of
dielectric elastomer actuators,” Applied Physics Letters, vol. 97, no. 16, p. 162908, 2010.
[123] B. Li, J. Zhou, and H. Chen, “Electromechanical stability in charge-controlled dielectric elastomer
actuation,” Applied Physics Letters, vol. 99, no. 24, p. 244101, 2011.
[124] Z. Alameh, S. Yang, Q. Deng, and P. Sharma, “Emergent magnetoelectricity in soft materials,
instability, and wireless energy harvesting,” Soft Matter, vol. 14, pp. 5856–5868, 2018.
[125] Q. Deng, L. Liu, and P. Sharma, “Flexoelectricity in soft materials and biological membranes,”
Journal of the Mechanics and Physics of Solids, vol. 62, pp. 209–227, 2014.
[126] Q. Deng, L. Liu, and P. Sharma, “Electrets in soft materials: Nonlinearity, size eﬀects, and giant
electromechanical coupling,” Physical Review E, vol. 90, no. 1, p. 012603, 2014.
[127] J. M. Ball, “Discontinuous equilibrium solutions and cavitation in nonlinear elasticity,” Phil. Trans.
R. Soc. Lond. A, vol. 306, no. 1496, pp. 557–611, 1982.
[128] C. Keplinger, T. Li, R. Baumgartner, Z. Suo, and S. Bauer, “Harnessing snap-through instability in
soft dielectrics to achieve giant voltage-triggered deformation,” Soft Matter, vol. 8, no. 2, pp. 285–
288, 2012.
[129] T. Li, C. Keplinger, R. Baumgartner, S. Bauer, W. Yang, and Z. Suo, “Giant voltage-induced
deformation in dielectric elastomers near the verge of snap-through instability,” Journal of the
Mechanics and Physics of Solids, vol. 61, no. 2, pp. 611–628, 2013.
114
Yang & Sharma: A tutorial on the stability and bifurcation Page 115
[130] A. Dorfmann and R. Ogden, “Nonlinear electroelastic deformations,” Journal of Elasticity, vol. 82,
no. 2, pp. 99–127, 2006.
[131] X. He, H. Yong, and Y. Zhou, “The characteristics and stability of a dielectric elastomer spherical
shell with a thick wall,” Smart Materials and Structures, vol. 20, no. 5, p. 055016, 2011.
[132] L. Dorfmann and R. W. Ogden, “Nonlinear response of an electroelastic spherical shell,” International
Journal of Engineering Science, vol. 85, pp. 163–174, 2014.
[133] Y.-X. Xie, J.-C. Liu, and Y. Fu, “Bifurcation of a dielectric elastomer balloon under pressurized
inﬂation and electric actuation,” International Journal of Solids and Structures, vol. 78, pp. 182–
188, 2016.
[134] F. Wang, C. Yuan, T. Lu, and T. Wang, “Anomalous bulging behaviors of a dielectric elastomer
balloon under internal pressure and electric actuation,” Journal of the Mechanics and Physics of
Solids, vol. 102, pp. 1–16, 2017.
[135] S. Rudykh, K. Bhattacharya, et al., “Snap-through actuation of thick-wall electroactive balloons,”
International Journal of Non-Linear Mechanics, vol. 47, no. 2, pp. 206–209, 2012.
[136] E. Bortot, “Analysis of multilayer electro-active spherical balloons,” Journal of the Mechanics and
Physics of Solids, vol. 101, pp. 250–267, 2017.
[137] W. T. Koiter, On the stability of elastic equilibrium (in Dutch with English summary). PhD thesis,
1945.
[138] J. Hutchinson and W. Koiter, “Postbuckling theory,” Appl. Mech. Rev, vol. 23, no. 12, pp. 1353–
1366, 1970.
[139] R. Casciaro, Computational asymptotic post-buckling analysis of slender elastic structures, pp. 195–
276. Vienna: Springer Vienna, 2005.
[140] W. Koiter, “Elastic stability, buckling and post-buckling behaviour,” in Proceedings of the IUTAM
Symposium on Finite Elasticity, pp. 13–24, Springer, 1981.
[141] G. J. Simitses, “Buckling and postbuckling of imperfect cylindrical shells: a review,” Applied
Mechanics Reviews, vol. 39, no. 10, pp. 1517–1524, 1986.
[142] N. Triantafyllidis and R. Peek, “On stability and the worst imperfection shape in solids with
nearly simultaneous eigenmodes,” International Journal of Solids and Structures, vol. 29, no. 18,
pp. 2281–2299, 1992.
[143] R. Peek and M. Kheyrkhahan, “Postbuckling behavior and imperfection sensitivity of elastic structures
by the lyapunov-schmidt-koiter approach,” Computer methods in applied mechanics and engineering,
vol. 108, no. 3-4, pp. 261–279, 1993.
115
Yang & Sharma: A tutorial on the stability and bifurcation Page 116
[144] G. Pampolini and N. Triantafyllidis, “Continuum electromechanical theory for nematic continua
with application to freedericksz instability,” Journal of Elasticity, vol. 132, pp. 219–242, 2018.
[145] Y. Cao and J. W. Hutchinson, “From wrinkles to creases in elastomers: the instability and
imperfection-sensitivity of wrinkling,” Proc. R. Soc. A, vol. 468, no. 2137, pp. 94–115, 2012.
[146] J. W. Hutchinson, “The role of nonlinear substrate elasticity in the wrinkling of thin ﬁlms,” Philosophical
Transactions of the Royal Society A: Mathematical, Physical and Engineering Sciences,
vol. 371, no. 1993, p. 20120422, 2013.
[147] C. Chiang Foo, S. Cai, S. Jin Adrian Koh, S. Bauer, and Z. Suo, “Model of dissipative dielectric
elastomers,” Journal of Applied Physics, vol. 111, no. 3, p. 034102, 2012.
[148] C. Chiang Foo, S. Jin Adrian Koh, C. Keplinger, R. Kaltseis, S. Bauer, and Z. Suo, “Performance of
dissipative dielectric elastomer generators,” Journal of Applied Physics, vol. 111, no. 9, p. 094107,
2012.
[149] S. Wang, M. Decker, D. L. Henann, and S. A. Chester, “Modeling of dielectric viscoelastomers with
application to electromechanical instabilities,” Journal of the Mechanics and Physics of Solids,
vol. 95, pp. 213 – 229, 2016.
[150] Y. Xiao and K. Bhattacharya, “A continuum theory of deformable, semiconducting ferroelectrics,”
Archive for Rational Mechanics and Analysis, vol. 189, pp. 59–95, 2008.
[151] F. Darbaniyan, K. Dayal, L. Liu, and P. Sharma, “Designing soft pyroelectric and electrocaloric
materials using electrets,” Soft Matter, vol. 15, pp. 262–277, 2019.
[152] E. Cerda, S. Chaieb, F. Melo, and L. Mahadevan, “Conical dislocations in crumpling,” Nature,
vol. 401, no. 6748, p. 46, 1999.
[153] E. Cerda and L. Mahadevan, “Conical surfaces and crescent singularities in crumpled sheets,”
Physical Review Letters, vol. 80, no. 11, p. 2358, 1998.
[154] P. Kodali, G. Saravanavel, and S. Sambandan, “Crumpling for energy: Modeling generated power
from the crumpling of polymer piezoelectric foils for wearable electronics,” Flexible and Printed
Electronics, vol. 2, no. 3, p. 035005, 2017.
[155] B. Wang, S. Yang, and P. Sharma, “Flexoelectricity as a universal mechanism for energy harvesting
from crumpling of thin sheets,” Physical Review B, vol. 100, p. 035438, 2019.
[156] P. P. Casta˜neda and M. Siboni, “A ﬁnite-strain constitutive theory for electro-active polymer
composites via homogenization,” International Journal of Non-Linear Mechanics, vol. 47, no. 2,
pp. 293 – 306, 2012. Nonlinear Continuum Theories.
[157] M. H. Siboni and P. P. Casta˜neda, “Fiber-constrained, dielectric-elastomer composites: Finitestrain
response and stability analysis,” Journal of the Mechanics and Physics of Solids, vol. 68,
pp. 211 – 238, 2014.
116
Yang & Sharma: A tutorial on the stability and bifurcation Page 117
[158] M. H. Siboni, R. Avazmohammadi, and P. P. Casta, “Electromechanical instabilities in ﬁberconstrained,
dielectric-elastomer composites subjected to all-around dead-loading,” Mathematics
and Mechanics of Solids, vol. 20, no. 6, pp. 729–759, 2014.
[159] M. H. Siboni and P. P. Casta˜neda, “Fiber-constrained dielectric elastomer composites: Finite
deformation response and instabilities under non-aligned loadings,” International Journal of Solids
and Structures, 2019.
[160] S. Rudykh and G. deBotton, “Stability of anisotropic electroactive polymers with application to
layered media,” Zeitschrift fr angewandte Mathematik und Physik, vol. 62, no. 6, pp. 1131–1142,
2011.
[161] S. Rudykh, K. Bhattacharya, and G. deBotton, “Multiscale instabilities in soft heterogeneous
dielectric elastomers,” Proceedings of the Royal Society A: Mathematical, Physical and Engineering
Sciences, vol. 470, no. 2162, p. 20130618, 2014.
[162] S. A. Spinelli and O. Lopez-Pamies, “Some simple explicit results for the elastic dielectric properties
and stability of layered composites,” International Journal of Engineering Science, vol. 88, pp. 15
– 28, 2015. Special Issue on.
[163] D. K. Vu, P. Steinmann, and G. Possart, “Numerical modelling of non-linear electroelasticity,”
International Journal for Numerical Methods in Engineering, vol. 70, no. 6, pp. 685–704, 2007.
[164] D. Vu and P. Steinmann, “On 3-d coupled bem–fem simulation of nonlinear electro-elastostatics,”
Computer Methods in Applied Mechanics and Engineering, vol. 201-204, pp. 82 – 90, 2012.
[165] C. Miehe, D. Rosato, and B. Kiefer, “Variational principles in dissipative electro-magnetomechanics:
A framework for the macro-modeling of functional materials,” International Journal
for Numerical Methods in Engineering, vol. 86, no. 10, pp. 1225–1276, 2011.
[166] D. Z¨ah and C. Miehe, “Computational homogenization in dissipative electro-mechanics of functional
materials,” Computer Methods in Applied Mechanics and Engineering, vol. 267, pp. 487 –
510, 2013.
[167] C. Miehe, D. Vallicotti, and S. Teichtmeister, “Homogenization and multiscale stability analysis in
ﬁnite magneto-electro-elasticity,” Gesellschaft f. Angewandte Mathematik und Mechanik, vol. 38,
pp. 313–343, 2015.
[168] C. Miehe, D. Vallicotti, and S. Teichtmeister, “Homogenization and multiscale stability analysis in
ﬁnite magneto-electro-elasticity. application to soft matter ee, me and mee composites,” Computer
Methods in Applied Mechanics and Engineering, vol. 300, pp. 294 – 346, 2016.
[169] E. Polukhov, D. Vallicotti, and M.-A. Keip, “Computational stability analysis of periodic electroactive
polymer composites across scales,” Computer Methods in Applied Mechanics and Engineering,
vol. 337, pp. 165 – 197, 2018.
117
Yang & Sharma: A tutorial on the stability and bifurcation Page 118
[170] M. Winterhalter and W. Helfrich, “Deformation of spherical vesicles by electric ﬁelds,” Journal of
Colloid and Interface Science, vol. 122, pp. 583–586, 1988.
[171] M. Kummrow and W. Helfrich, “Deformation of giant lipid vesicles by electric ﬁelds,” Physical
Review A, vol. 44, pp. 8356–8360, 1991.
[172] T. Yamamoto, S. Aranda-Espinoza, R. Dimova, and R. Lipowsky, “Stability of spherical vesicles
in electric ﬁelds,” Langmuir, vol. 26, pp. 12390–12407, 2010.
[173] T. Portet, C. Mauroy, V. D´emery, T. Houles, J.-M. Escoﬀre, D. S.Dean, and M.-P. Rols, “Destabilizing
giant vesicles with electric ﬁelds: An overview of current applications,” The Journal of
Membrane Biology, vol. 245, pp. 555–564, 2012.
[174] P. M. Vlahovska, “Electrohydrodynamics of drops and vesicles,” Annual Review of Fluid Mechanics,
vol. 51, no. 1, pp. 305–330, 2019.
118