# Universal deformations of ideal liquid crystal elastomers

Victoria Lee  Saint-Gobain Competency Research Laboratory, Northborough MA
01532 222This work was conducted while Victoria Lee was affiliated with the
California Institute of Technology. Kaushik Bhattacharya111Corresponding
Author. Email<EMAIL_ADDRESS>California Institute of Technology, Pasadena
CA 91125

###### Abstract

Liquid crystal elastomers are rubber-like solids with liquid crystalline
mesogens (stiff, rod-like molecules) incorporated either into the main chain
or as a side chain of the polymer. These solids display a range of unusual
thermo-mechanical properties as a result of the coupling between the entropic
elasticity of rubber and the orientational phase transitions of liquid
crystals. One of these intriguing properties is the soft behavior, where it is
able to undergo significant deformations with almost no stress. While the
phenomenon is well-known, it has largely been examined in the context of
homogenous deformations. This paper investigates soft behavior in complex
inhomogeneous deformations. We model these materials as hyperelastic,
isotropic, incompressible solids and exploit the seminal work of Ericksen, who
established the existence of non-trivial universal deformations, those that
satisfy the equations of equilibrium in every hyperelastic, isotropic,
incompressible solid. We study the inflation of spherical and cylindrical
balloons, cavitation and bending.

Dedicated to the memory of Jerald L. Ericksen, an original thinker and
inspiring teacher

##  1 Introduction

Liquid crystal elastomers are rubber-like solids with liquid crystalline
mesogens (stiff, rod-like molecules) incorporated either into the main chain
or as a side chain of the polymer. These solids combine the entropic
elasticity of rubber with the orientational phase transitions of liquid
crystals. The deformation and orientational order are coupled, giving rise to
a whole host of very interesting thermo-mechanical properties. Originally
envisioned by de Gennes in 1975 [13], and after early attempts at synthesis
[17], these materials were first reliably synthesized by Küpfer and Finkelmann
in 1991 [23]. The development of a statistical mechanical theory [8], an
understanding of the mathematical structure of this theory [15], the
recognition that the genesis, the state in which the polymer is cross-linked,
plays a critical role [7, 36], the discovery of new chemistries [39, 40] and a
variety of directed methods of synthesis [37, 3] have made these materials
widely available, and the subject of both fundamental and applied studies. We
refer the reader to Warner and Terentjev [38] for a comprehensive
introduction, and White and Broer [39] for a recent review on thermo-
mechanical applications.

The limited cross-linking of liquid crystal elastomers enables the liquid
crystalline mesogens to undergo their usual order-disorder phase transition
with changes in temperature. They are disordered (randomly oriented) at high
temperatures but develop nematic order where they are locally oriented in a
particular direction (represented by a unit vector or director) at lower
temperature due to steric interactions. However, the coupling to the elastomer
results in a change of shape, an elongation along the nematic director and
transverse contraction, as the liquid crystal elastomer is cooled from its
high-temperature isotropic state to its low-temperature nematic state. This
change of shape can be exploited for actuation. Further, mechanical loads can
reorient the director in some specimens, and this leads to an intriguing soft
behavior. This soft behavior is the focus of the current paper.

The soft behavior was first observed by Finkelmann and coworkers [23, 22] in a
carefully synthesized liquid crystal elastomer sheet with uniform initial
nematic director. They subjected these sheets to uniaxial tension
perpendicular to the initial nematic and observed that the the sheet could
undergo stretches as large as 3 (up to 200% strain) with little stress. In-
situ observations using birefringence microscopy established that the soft
behavior is related to the reorientation of the director. Intriguingly, the
reorientation is not uniform: instead the sheet is divided into ‘stripe
domains’, stripes on the scale of microns, with the director rotating
clockwise and counter-clockwise in alternate stripes. Bladon, Warner and
Terentjev [8] developed the so-called neo-classical theory of liquid crystal
elastomers and used it to explain the stripe domains and soft-behavior.
DeSimone and Dolzmann [14] showed that the Bladon, Warner and Terentjev energy
was not (rank-one) convex, and the stripe domains were a natural consequence
of this non-convexity. Subsequently, they computed the full relaxation of this
energy and showed the possibility of a rich class of soft deformations [15].
More recently, the soft behavior has been investigated in biaxial stretch
revealing a remarkable liquid-like in-plane behavior [9, 34, 42]. This soft
behavior is the basis of a number of proposed applications including impact
resistance [32]. Still, there is limited study of soft behavior under
inhomogeneous deformations.

An isotropic-genesis material is isotropic, and thus the symmetry-breaking
isotropic-to-nematic phase transition leads to a (infinite number of) nematic
states or variants. In particular, all director orientations are equivalent,
and thus the material may undergo spontaneous stretch in any direction.
Further, they are free to reorient to accommodate any imposed deformation.
Finally, different regions can have different directors, leading to the
formation of domains. However, the domains cannot be arbitrary: they have to
satisfy mechanical (Hadamard) compatibility conditions across the domain walls
(boundaries across which directors suffer a jump). This compatibility
condition gives rise to the stripe domains, but also allows other more complex
patterns [14, 15]. In practice, there is some disorder in the cross-link
density that adds some local random anisotropy [27]. This is often described
as ‘non-ideality’ and leads to what is described as ‘semi-soft’ behavior,
where one needs a small stress to reorient the director [11, 10, 6]. Still
neglecting this non-ideality provides significant amounts of insight into the
behavior of a liquid crystal elastomer.

In this paper, we study how soft elasticity of liquid crystal elastomers
manifests itself during complex inhomogeneous deformations and affects the
response of structures. We do so by exploiting the seminal work of Ericksen
[16]. He showed that there exist certain (parameterized) families of
inhomogeneous deformations that automatically satisfy the equation of
mechanical equilibrium in every isotropic, incompressible, hyperelastic body.
Briefly, incompressibility leads to a hydrostatic pressure that is not
constitutively determined, and this indeterminacy in the pressure enables the
satisfaction of the equilibrium equation. Remarkably, the families of
deformation are incredibly rich and encompass many common situations
encountered in application. Since these deformations automatically satisfy the
equation of equilibrium, one only has to determine the parameters, and one can
do so from macroscopic equilibrium and boundary conditions. This remarkable
result explains the success of the semi-inverse method of Rivlin [29, 28, 30,
31] and is the foundation of much work in finite elasticity since (see [24,
43] for recent applications, and [41] for a discussion and extension to
inhomogeneous bodies).

An isotropic-genesis liquid crystal elastomer in the ideal limit (described by
the neo-classical theory of Bladon, Warner and Terentjev [8] and its
relaxation by DeSimone and Dolzmann [15]) is an isotropic, incompressible
hyperelastic solid. Therefore, the universal deformations of Ericksen satisfy
the equations of equilibrium. The fact that it is soft is irrelevant, and
therefore these universal deformations are an ideal avenue to probe the soft
elasticity of liquid crystal elastomers in inhomogeneous deformations.

We study three families of inhomogeneous deformation. The first family is
spherically symmetric expansion/compression, and we study two problems. The
first problem is the inflation of a (possibly thick-walled) spherical balloon.
Spherical balloons subjected to internal pressure undergo the “balloon
instability”, where the radius changes in a discontinuous manner with
increasing pressure [2]. We show that the soft behavior reduces the critical
pressure at which this instability occurs and vastly enhances the radius jump.
This problem has also been recently studied by Giudici and Biggins [19].

Figure 1: Schematic showing the cross-section of a disk of nematic elastomer
bonded to parallel plates, which are stretched in uniaxial tension.

The second problem in the spherical family is cavitation, motivated by the
pioneering experiments and analysis of Gent and Lindley [18]. They study a
short, rubber cylinder that is glued to grips and subjected to tension, as
shown schematically in Figure 1. They observe an unexpected instability at a
critical load accompanied by internal ruptures. The interior of the cylinder
is subjected to uniform hydrostatic tension before rupture. So they study a
spherical annulus subjected to external hydrostatic tension; they find by
using a neo-Hookean material model that the inner void expands uncontrollably
at a critical tension, the critical tension is independent of the initial
radius of the inner void, and it corresponds to the stress associated with the
instability. They suggest that any defect would grow rapidly at this critical
tension and eventually lead to the nucleation of an internal rupture. Ball [5]
showed that an initial hole is not necessary, and discontinuous solutions can
arise as a result of a lack of growth in the energy density. We follow Gent
and Lindley [18] and study the expansion of a spherical annulus subjected to
external hydrostatic tension. We show that the soft behavior promotes
cavitation. Mihai and Goriely [26] studied this problem but from the point of
view of uncertainty, and how the instability is affected by the anisotropy
parameter and shear modulus.

The second family of deformations we study is cylindrically symmetric
inflation/compression, and apply it to the problem of inflation of a
cylindrical balloon. We again have the balloon instability, and the soft
behavior reduces the critical pressure at which this instability occurs and
vastly enhances the radius jump. A closely related problem, instability of a
pressurized thin-walled cylindrical balloon under axial loads, has been
studied by He et al. [20]. Further, the problem of inflation of a nematic-
genesis cylindrical balloon and its potential use as a pump is presented
elsewhere [25]. The third family concerns bending. We show how the soft
behavior manifests itself as a plateau in the moment-curvature relationship.
In each of these families, we show how the underling domain patterns would
evolve. A fourth family involving the torsion of a cylinder has been studied
by Baardink and Cesana [4].

The paper is organized as follows. We begin by presenting the material model
in Section 2. The Bladon-Terentjev-Warner [8] theory, like the neo-Hookean
constitutive model, is based on a Gaussian approximation to the statistics of
polymer chains, and as such is not appropriate for large deformations.
Therefore, we propose a generalization based on a generalized Mooney-Rivlin
constitutive model. This is similar to a model introduced by Agostiniani-
DeSimone [1]. We then present its relaxation and compute the stress associated
with this relaxed energy. We introduce Ericksen’s universal deformations in
Section 3. The following three sections study the three families of
deformation described above.

##  2 Constitutive relations

###  2.1 Stored energy density

Consider a specimen of the isotropic-genesis liquid crystal elastomer in the
stress-free isotropic state as the reference configuration that occupies the
domain $\Omega$. The material in the current state is nematic with a director
$\bm{n}$, $|\bm{n}|=1$, and deformation gradient $\bm{F}$ relative to the
reference state. Following Bladon-Terentjev-Warner, we postulate that the
material is incompressible ($\det\bm{F}=1$), and the stored energy density is

$W(\bm{F},\bm{n})=W^{\text{el}}(\bm{\ell}_{n}^{-1/2}\bm{F}),$ (1)

where

$\bm{\ell}_{n}=r^{-1/3}\left(\textbf{I}+(r-1){\bm{n}}\otimes{\bm{n}}\right)$
(2)

is the step-length tensor, and $r>1$ is a (temperature-dependent) parameter
that describes the degree of nematic order. Note that $\det\ \bm{\ell}_{n}=1$.
It is easy to verify that the energy is frame-indifferent under a change of
frame $\bm{x}\mapsto\bm{Q}\bm{x}+\bm{c}$, since
$\bm{n}\mapsto\bm{Q}\bm{n},\bm{F}\mapsto\bm{Q}\bm{F}$ leaves $W$ invariant for
frame-indifferent $W^{\text{el}}$. Similarly, since
$\bm{n}\mapsto\bm{n},\bm{F}\mapsto\bm{F}\bm{R}$ under a change of material
frame $\bm{X}\mapsto\bm{R}\bm{X}$, $W$ is isotropic for isotropic
$W^{\text{el}}$. We take $W^{\text{el}}$ to be the generalized Mooney-Rivlin
energy; so

$W(\bm{F},\bm{n})=\sum_{i=1}^{n}c_{i}(\tilde{I}_{1}-3)^{p_{i}}+\sum_{j=1}^{m}d_{j}(\tilde{I}_{2}-3)^{q_{j}},$
(3)

where

$\tilde{I}_{1}=\text{tr }\tilde{\bm{b}},\quad\tilde{I}_{2}=\text{tr }\text{cof
}\tilde{\bm{b}}=\text{tr
}(\tilde{\bm{b}})^{-1},\quad\tilde{\bm{b}}=\bm{\ell}_{n}^{-1/2}\bm{F}\bm{F}^{T}\bm{\ell}_{n}^{-1/2}.$
(4)

We take

$c_{i}\geq 0,\ p_{i}\geq 1,\ d_{j}\geq 0,\ q_{j}\geq 1,\quad i=1,\dots
n,j=1,\dots,m$ (5)

so that the corresponding $W^{\text{el}}$ is polyconvex. Note that the BTW
theory corresponds to $n=1,p_{1}=1,m=0$.

Agostiniani and DeSimone [1] introduced a slightly different generalization,

$W^{\text{AD}}(\bm{F},\bm{n})=\sum_{i=1}^{n}c_{i}(\tilde{I}_{1}^{p_{i}}-3)+\sum_{j=1}^{m}d_{j}(\tilde{I}_{2}^{q_{j}}-3).$
(6)

It is common to minimize $\bm{n}$ out and define

$\widehat{W}(\bm{F})=\min_{|\bm{n}|=1}W(\bm{F},\bm{n}).$ (7)

It is a long but straightforward calculation [15, 1] to show that

$\widehat{W}(\bm{F})=\widetilde{W}(s,t)=\begin{cases}\begin{aligned}
&\sum_{i=1}^{n}c_{i}\bigg{\lvert}r^{1/3}\left(\frac{s^{2}}{r}+\frac{t^{2}}{s^{2}}+\frac{1}{t^{2}}\right)-3\bigg{\rvert}^{p_{i}}\\\
&\quad+\sum_{j=1}^{m}d_{j}\bigg{\lvert}r^{-1/3}\left(\frac{r}{s^{2}}+\frac{s^{2}}{t^{2}}+t^{2}\right)-3\bigg{\rvert}^{q_{j}}\end{aligned}&\
\det\bm{F}=1\\\ \infty&\ \text{else}\end{cases},$ (8)

where $s$ is the largest singular value of $\bm{F}$, and $t$ is the largest
singular value of $\text{cof }\bm{F}$. In other words, if
$\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}$ are the ordered singular values of
$\bm{F}$, then $s=\lambda_{1}$ and $t=\lambda_{1}\lambda_{2}$. Unfortunately,
$\widehat{W}$ is not rank-one convex, and hence not quasiconvex. We refer the
reader to Dacorogna [12] for a broad introduction to convexity conditions and
their consequences for the existence of energy minimizers.

It is easy to follow the arguments of DeSimone and Dolzmann [15] as well as
Agostiniani and DeSimone [1] to compute the relaxation or quasiconvexification
of the $\widehat{W}$ to be

${\overline{W}}({\mathbf{F}})=W^{qc}(s,t)=\sum_{i=1}^{M}c_{i}\left(W_{1}(s,t)\right)^{p_{i}}+\sum_{j=1}^{N}d_{j}\left(W_{2}(s,t)\right)^{q_{j}},$
(9)

where

$W_{1}(s,t)=\begin{cases}0&(s,t)\in L\\\
\frac{r^{1/3}}{t^{2}}+\frac{2t}{r^{1/6}}-3&(s,t)\in M\\\
r^{1/3}\left(\frac{s^{2}}{r}+\frac{t^{2}}{s^{2}}+\frac{1}{t^{2}}\right)-3&(s,t)\in
S\\\ \infty&\text{else}\end{cases}\ $ (10)

and

$W_{2}(s,t)=\begin{cases}0&(s,t)\in L\\\
r^{-1/3}t^{2}+\frac{2r^{1/6}}{t}-3&(s,t)\in M\\\
r^{-1/3}\left(\frac{r}{s^{2}}+\frac{s^{2}}{t^{2}}+t^{2}\right)-3&(s,t)\in S\\\
\infty&\text{else}\end{cases},$ (11)

and the regions $L$, $M$, and $S$ are given by

$\displaystyle L$ $\displaystyle=\\{(s,t):t\leq s^{2},\ t\geq\sqrt{s},\ t\leq
r^{1/6}\\},$ (12) $\displaystyle M$ $\displaystyle=\\{(s,t):t\geq r^{1/6},\
t\leq s^{2},\ t\geq r^{-1/2}s^{2}\\},$ $\displaystyle S$
$\displaystyle=\\{(s,t):t\geq\sqrt{s},\ t\leq r^{-1/2}s^{2}\\}.$

These regions are shown in Figure 2. Henceforth, we refer to (9) as the
relaxed generalized Mooney-Rivlin (RGMR) model and the special case where
$M=1,N=0$ as the relaxed Bladon-Terentjev-Warner (RBTW) model333We also have
$p_{1}=1$ in the original Bladon-Terentjev-Warner model but we do not require
it here..

Figure 2: Regions of $L$, $M$, and $S$ in the phase diagram of $(s,t)$.

We omit the proof of the relaxation since it closely follows that of
Agostiniani and DeSimone [1]. Instead we provide a brief overview. Recall that
$W$ is isotropic. It follows that $\widehat{W}$ obtained by minimizing over
all possible orientations and $\widetilde{W}$ obtained by relaxation are also
isotropic. Therefore, the energy can only depend on the principal values
$\lambda_{1}\geq\lambda_{2}\geq\lambda_{3}>0$ of the deformation gradient
$\bm{F}$. Incompressibility dictates that
$\lambda_{1}\lambda_{2}\lambda_{3}=1$, and thus the energies only depend on
$\lambda_{1},\lambda_{2}$. Since $s=\lambda_{1}$ and
$t=\lambda_{1}\lambda_{2}$, we may write the energies only in terms of $s,t$
as we do above. $\widetilde{W}=0$ if and only if $s=r^{1/3},t=r^{1/6}$, i.e.,
when $\lambda_{1}=r^{1/3},\lambda_{2}=\lambda_{3}=r^{-1/6}$ or when the
specimen is elongated in one directlon and laterally compressed. Further,
$\widetilde{W}$ grows away from this state.

Now, observe that $s$ and $t$ are polyconvex and non-negative ($s$ is convex
in $\bm{F}$ and $t$ in $\text{cof }\bm{F}$). Therefore, $\overline{W}$ is
polyconvex if it is non-decreasing in $s$ and $t$. Therefore, replacing the
decreasing parts of $\widetilde{W}$ with constants gives an upper bound on the
relaxation. However, this bound can be attained by laminates, and we obtain
the relaxation. In the region marked $L$, we relax in both $s,t$ and the bound
is attained by double laminates. Here the relaxed energy is identically zero.
In the region marked $M$, we relax only in $s$ and therefore the energy is
independent of $s$. The bound here is attained by laminates (with a unique
Young measure [9]). Finally, there is no relaxation in the region marked $S$.

Note that $\bm{F}=\bm{I}$ or $s=t=1$ belongs to $L$, and thus the
(macroscopically) undeformed state relative to the high-temperature state is
unstressed. As we deform, it is initially stress-free as $\bm{F}$ traverses
$L$, and it becomes stressed as $\bm{F}$ reaches $M$ and $S$ (two equal
principal Cauchy stress in $M$).

###  2.2 Stress

We can readily compute the Cauchy stress from the energy density above as

$\bm{\sigma}=-p\bm{I}+\sum_{i=1}^{3}\hat{\sigma}_{i}(s,t)\hat{\bm{v}}_{i}\otimes\hat{\bm{v}}_{i},\quad\text{where}\quad\hat{\sigma}_{i}(s,t)=\lambda_{i}\frac{\partial{W}}{\partial{\lambda_{i}}},$
(13)

$p$ is an unknown hydrostatic pressure to be determined from equilibrium
resulting from incompressibility, $\lambda_{i}$ are the ordered principal
stretches as above and $\hat{\bm{v}}_{i}$ are the principal directions
(eigenvalues and eigenvectors of the left Cauchy-Green tensor, or
$\bm{b}=\bm{F}{\bm{F}}^{T}=\sum_{i=1}^{3}\lambda_{i}^{2}\hat{\bm{v}}_{i}\otimes\hat{\bm{v}}_{i}$).
It is a long but straightforward calculation to conclude:

$\displaystyle\hat{\sigma}_{1}(s,t)$ $\displaystyle=$
$\displaystyle\begin{cases}0&(s,t)\in L\\\
\sum_{i=1}^{M}c_{i}p_{i}|W_{1}|^{p_{i}-1}2r^{-1/6}t+\sum_{j=1}^{N}d_{j}q_{j}|W_{2}|^{q_{j}-1}2r^{-1/3}t^{2}&(s,t)\in
M\\\
\sum_{i=1}^{M}c_{i}p_{i}|W_{1}|^{p_{i}-1}2r^{-2/3}s^{2}+\sum_{j=1}^{N}d_{j}q_{j}|W_{2}|^{q_{j}-1}2r^{-1/3}(t^{2}+s^{2}t^{-2})&(s,t)\in
S\end{cases},$ $\displaystyle\hat{\sigma}_{2}(s,t)$ $\displaystyle=$
$\displaystyle\begin{cases}0&(s,t)\in L\\\
\sum_{i=1}^{M}c_{i}p_{i}|W_{1}|^{p_{i}-1}2r^{-1/6}t+\sum_{j=1}^{N}d_{j}q_{j}|W_{2}|^{q_{j}-1}2r^{-1/3}t^{2}&(s,t)\in
M\\\
\sum_{i=1}^{M}c_{i}p_{i}|W_{1}|^{p_{i}-1}2r^{1/3}t^{2}s^{-2}+\sum_{j=1}^{N}d_{j}q_{j}|W_{2}|^{q_{j}-1}2(r^{-1/3}t^{2}+r^{2/3}s^{-2})&(s,t)\in
S\\\ \end{cases},$ (14) $\displaystyle\hat{\sigma}_{3}(s,t)$ $\displaystyle=$
$\displaystyle\begin{cases}0&(s,t)\in L\\\
\sum_{i=1}^{M}c_{i}p_{i}|W_{1}|^{p_{i}-1}2r^{1/3}t^{-2}+\sum_{j=1}^{N}d_{j}q_{j}|W_{2}|^{q_{j}-1}2r^{1/6}t^{-1}&(s,t)\in
M\\\
\sum_{i=1}^{M}c_{i}p_{i}|W_{1}|^{p_{i}-1}2r^{1/3}t^{-2}+\sum_{j=1}^{N}d_{j}q_{j}|W_{2}|^{q_{j}-1}2(r^{-1/3}s^{2}t^{-2}+r^{2/3}s^{-2})&(s,t)\in
S\end{cases},$

where $W_{1},W_{2}$ are as defined in (10, 11).

Note that the stress is purely hydrostatic in the region $L$ and has two equal
principal values $\sigma_{1}=\sigma_{2}$ in region $M$.

###  2.3 Homogeneous deformations

We desribe certain homogeneous deformations to understand the nature of the
constitutive relations. These are chosen to be representative of the states of
stretch and stress that we encounter in our examples later.

##### Uniaxial Stress

Here a long cylindrical specimen is subject to tensile stress along the axis
and the lateral surfaces are traction free. These experiments go back to
Finkelmann and coworkers [23, 22] in specially prepared monodomain specimens,
and to Urayama [36] in isotropic genesis polydomain specimens. Picking a
coordinate system so that the $1-$direction is along the axes, we have
$\sigma_{11}=\sigma,\sigma_{22}=\sigma_{33}=0$. We make the ansatz

$\bm{F}=\begin{pmatrix}\lambda&&\\\ &\lambda^{-1/2}&\\\
&&\lambda^{-1/2}\end{pmatrix}\implies\bm{\sigma}=\begin{pmatrix}\hat{\sigma}_{1}(\lambda,\sqrt{\lambda})-\hat{\sigma}_{3}(\lambda,\sqrt{\lambda})&&\\\
&0&\\\ &&0\end{pmatrix}$ (15)

where we have used the condition $\sigma_{33}=0$ and the fact that
$s=\lambda,t=\lambda^{1/2}$. We start at $\\{1,1\\}$ and traverse along the
lower boundary $t=\sqrt{s}$ of Figure 2. Therefore we start in region $L$ and
transition to the region $S$ when $\lambda=r^{1/3}$: we do not encounter the
region $M$. We substitute the values of $s,t$ in (14) to find the state of
stress. We plot this stress vs. stretch relation for both the RGMR and RBTW
models for the parameters $r=8$ and

$M=2,\ N=1,\ c_{1}=1.0\times 10^{5}\text{Pa},\ c_{2}=1.90\times
10^{2}\text{Pa},\ d_{1}=1.59\times 10^{-2}\text{Pa},\ p_{1}=1.3,\ p_{2}=2,\
q_{1}=2.$ (16)

in Figure 3(a). We have a soft behavior with material in the region $L$ with
microstructure and zero stress till $\lambda=r^{1/3}=2$; subsequently the
stress increases as the material enters region $S$ and there is no
microstructure. These are consistent with experimental observations [36]444See
[34, 42] for detailed studies using a model with non-ideality. . We do not see
a difference between the RGMR and RBTW models for the values of the stretch
that are plotted.

Figure 3: Stress-stretch behavior in homogeneous deformations (a) Uniaxial
tension. (b) Equibiaxial stretch. (c) Uniaxial Stretch. The dark black line,
blue continuous line and red dashed line indicates that the material is in
region L, M and S respectively. The dark line indicates the RGMR model while
the light line indicates the RBTW model.

##### Equibiaxial Stretch

Here a sheet with traction-free surfaces is stretched equally along two
perpendicular in-plane directions. This was studied experimentally by
Tokumuoto et al. [34]. Picking a coordinate system so that the $3-$direction
is normal to the sheet, we have
$\sigma_{33}=0,\lambda_{1}=\lambda_{2}=\lambda$. Therefore,

$\bm{F}=\begin{pmatrix}\lambda&&\\\ &\lambda&\\\
&&\lambda^{-2}\end{pmatrix}\implies\bm{\sigma}=\begin{pmatrix}\hat{\sigma}_{1}(\lambda,\lambda^{2})-\hat{\sigma}_{3}(\lambda,\lambda^{2})&&\\\
&\hat{\sigma}_{1}(\lambda,\lambda^{2})-\hat{\sigma}_{3}(\lambda,\lambda^{2})&\\\
&&0\end{pmatrix}$ (17)

where we have used the condition $\sigma_{33}=0$, and the fact that
$\hat{\sigma}_{1}=\hat{\sigma}_{2}$ and $s=\lambda,t=\lambda^{2}$. We again
start at $\\{1,1,\\}$ and traverse along the upper boundary $t=s^{2}$ of
Figure 2. Therefore we start in region $L$ and transition to the region $M$
when $\lambda=r^{1/12}$: we do not encounter the region $S$. We substitute the
values of $s,t$ in (14) to find the state of stress. We plot this for both the
RGMR and RBTW models for the parameters $r=8$ and (16) and this shown in
Figure 3(b). We have a soft behavior with zero stress till
$\lambda=r^{1/12}\approx 1.18$ and then it increases as the material enters
region $M$. We do not enter the region $S$ and the material always has
microstructure. These are consistent with experimental observations
[34]${}^{\ref{foot:nonideal}}$. We do not see a difference between the RGMR
and RBTW models for small values of stretch, but they eventually diverge for
larger values of stretch. These become relevant in our studies of balloons
later.

##### Uniaxial Stretch

Here a sheet with traction-free surfaces is stretched along one in-plane
direction while being constrained in the other. This was studied
experimentally by Tomuoto et al. [34]. Picking a coordinate system so that the
$3-$direction is normal to the sheet, we have
$\sigma_{33}=0,\lambda{1}=\lambda,=\lambda{2}=1$. Therefore,

$\bm{F}=\begin{pmatrix}\lambda&&\\\ &1&\\\
&&\lambda^{-1}\end{pmatrix}\implies\bm{\sigma}=\begin{pmatrix}\hat{\sigma}_{1}(\lambda,\lambda)-\hat{\sigma}_{3}(\lambda,\lambda)&&\\\
&\hat{\sigma}_{2}(\lambda,\lambda)-\hat{\sigma}_{3}(\lambda,\lambda)&\\\
&&0\end{pmatrix}.$ (18)

where we have used the condition $\sigma_{33}=0$ and the fact that
$s=t=\lambda$. We start at $\\{1,1\\}$ and traverse along the diagonal of
Figure 2. Therefore we start in region $L$, transition to the region $M$ when
$\lambda=r^{1/6}$ and again transition to the region $S$ when
$\lambda=r^{1/2}$. We substitute the values of $s,t$ in (14) to find the state
of stress. We plot this for both the RGMR and RBTW models for the parameters
$r=8$ and (16) and this shown in Figure 3(c). We have a soft behavior with
zero stress till $\lambda=r^{1/6}=\sqrt{2}$ and then it increases as the
material enters region $M$ and the stress begins to rise. However, remarkably,
$\sigma_{1}=\sigma_{2}$ even though $\lambda_{1}\neq\lambda_{2}$. In other
words, we have shear strain but no shear stress, and therefore this has been
described as “in-plane liquid like behavior” [34]. We transition out of region
$M$ to region $S$ at $\lambda=r^{1/2}\approx 2.83$ and the two components of
stress diverge. These are consistent with experimental observations
[34]${}^{\ref{foot:nonideal}}$. We do not see a difference between the RGMR
and RBTW models for small values of stretch, but they eventually diverge for
larger values of stretch.

##  3 Ericksen’s universal deformations

In a seminal work, Ericksen [16] showed that it is possible to find rich
classes of non-trivial, inhomogeneous555Homogeneous deformations satisfy the
equilibrium equation trivially. solutions to the equilibrium equations that
hold for every isotropic, incompressible, hyperelastic solid (i.e.,
independent of the specific constitutive relation) in the absence of body
forces. Ericksen identified four families of inhomogeneous universal
deformations. A fifth family was independently discovered by Kingbell and
Shield [21] and Singh and Pipkin [33].

It is convenient to use different coordinate systems to describe deformation;
we use upper-case letters to describe the reference configuration ($(X,Y,Z)$,
$(R,\Theta,Z)$ or $(R,\Theta,\Phi)$ for rectangular Cartesian, cylindrical or
spherical respectively) and lower-case letters to describe the current
configuration ($(x,y,z)$, $(\rho,\theta,z)$ or $(\rho,\theta,\phi)$).
Constants $a,b,c,d,e,f$ parameterize the solutions. The five families are:
Family 1: Bending, stretching and shearing of a rectangular block

$\rho=\sqrt{2aX+d},\theta=bY,z=\frac{Z}{ab}-bcY.$ (19)

Family 2: Straightening, stretching and shearing of a sector of a tube

$x=\frac{1}{2}ab^{2}R^{2},y=\frac{\Theta}{ab},z=\frac{Z}{b}-\frac{c\Theta}{ab}.$
(20)

Family 3: Inflation, bending, torsion, extension and shearing of an annular
wedge, with $a(cf-de)=1$

$\rho=\sqrt{aR^{2}+b},\theta=c\Theta+dZ,z=e\Theta+fZ.$ (21)

Family 4: Inflation or eversion of a sector of a spherical shell

$\rho=[sR^{3}+a]^{1/3},\theta=s\Theta,\phi=\Phi,\quad s=\pm 1.$ (22)

Family 5: Inflation, bending, extension and azimuthal shearing of an annular
wedge [21, 33]

$\rho=a^{1/2}R,\theta=d\ln{(bR)}+c\Theta,z=eZ,\ ace=1.$ (23)

It turns out that these families are also almost exhaustive: one interesting
mathematical possibility remains open, but it is unknown if that possibility
has any non-trivial deformations.

We note that while these universal deformations allow us to solve boundary-
value problems as described above, they are not guaranteed to be the unique
solution. However, experience suggests that these are often unique for small
distortion, but other solutions bifurcate from this branch at larger
distortion.

Since ideal LCEs are isotropic, incompressible, hyperelastic solids, we can
exploit these universal deformations to study some boundary-value problems. We
focus on three families – 1, 3 and 4.

##  4 Radial deformation of a spherical shell

###  4.1 Radial deformation of a spherical shell

Figure 4: (a) A spherical shell subjected to internal and external pressure.
(b) Three different cases in a spherical nematic shell (see (30) for
$\lambda_{i,o}^{\pm})$.

####  4.1.1 Kinematics

We consider a spherical shell subjected to an internal and external pressure
as shown in Figure 4(a), and consider a special case of Family 4 with $s=+1$:

$\rho=(R^{3}+a)^{1/3},\theta=\Theta,\phi=\Phi.$ (24)

It is easy to verify that the deformation gradient and left Cauchy-Green
tensor are

$\bm{F}=\begin{pmatrix}\frac{1}{\lambda^{2}}&&\\\ &\lambda&\\\
&&\lambda\end{pmatrix},\quad\bm{b}=\begin{pmatrix}\frac{1}{\lambda^{4}}&&\\\
&\lambda^{2}&\\\
&&\lambda^{2}\end{pmatrix},\quad\text{where}\quad\lambda={(R^{3}+a)^{1/3}\over
R}$ (25)

in spherical coordinates with $a>0$ for inflation. It follows that the
principal stretches are $\lambda,\lambda,1/\lambda^{2}$ with the principal
stretches in the $\theta$, $\phi$ and $\rho$ directions respectively. Further,
$s=\lambda,t=\lambda^{2}$. Since $t=s^{2}$, each element is subjected to equi-
biaxial stretch normal to the radius. We only expect to encounter the regions
$L$ and $M$ – see Figure 2 with the transition when $\lambda=r^{1/12}$.

Further, $\lambda$ is monotonically decreasing in $R$, and therefore we expect
region $L$ in the interior and $M$ in the exterior, with the transition taking
place at the radius $R^{*}$, where $\lambda(R^{*})=r^{1/12}$:

$L=\\{R\geq R^{*}\\},\quad M=\\{R<R^{*}\\},\quad\text{where}\quad
R^{*}=a^{1/3}(r^{1/4}-1)^{-1/3}.$ (26)

Depending on the constant $a$ and the inner and outer radii, we have three
possible cases:

  * •

Case L: The entire shell is in region $L$ when

$\lambda_{i}\leq r^{1/12}\iff\lambda_{i,o}\leq\lambda_{i,o}^{-}\iff R^{*}\leq
R_{i};$ (27)

  * •

Case ML: The inner region of the shell $R\in[R_{i},R^{*}]$ is in $M$, and the
outer region $R\in[R^{*},R_{o}]$ is in $L$ when

$\lambda_{o}<r^{1/12}<\lambda_{i}\iff\lambda_{i,o}^{-}<\lambda_{o}<\lambda_{i,o}^{+}\iff
R_{i}<R^{*}<R_{o};$ (28)

  * •

Case M: The entire shell is in region $M$ when

$r^{1/12}\leq\lambda_{o}\iff\lambda_{i,o}^{+}\leq\lambda_{i,o}\iff R_{o}\leq
R^{*}.$ (29)

Above, $\lambda_{i}=\lambda(R_{i}),\lambda_{o}=\lambda(R_{o})$, and

$\lambda_{i}^{-}=r^{1/12},\quad\lambda_{i}^{+}={(R_{i}^{3}+R_{o}^{3}(r^{1/4}-1))^{1/3}\over
R_{i}},\quad\lambda_{o}^{-}={(R_{o}^{3}+R_{i}^{3}(r^{1/4}-1))^{1/3}\over
R_{o}},\quad\lambda_{o}^{+}=r^{1/12}.$ (30)

A diagram illustrating the various cases can be seen in Figure 4(b), and a
table of the values $\lambda_{0}^{\pm}$ for various $t$ is given in Table 1.

$r$ | $\lambda_{0}^{-}$ | $\lambda_{0}^{+}$
---|---|---
2 | 1.045 | 1.059
4 | 1.095 | 1.122
6 | 1.125 | 1.161
8 | 1.148 | 1.189
Table 1: Transitions between different regimes

In summary, an ideal nematic shell always has microstructure with all parts of
the shell in either region $L$ or $M$; when there are two regions, the inner
layer is in the liquid-like region, and the outer layer is in the
microstructure region.

####  4.1.2 Equilibrium

The equations of equilibrium in spherical coordinates are

$\begin{aligned}
\rho&:\frac{\partial\sigma_{\rho\rho}}{\partial\rho}+2{\sigma_{\rho\rho}\over\rho}+{1\over\rho}\frac{\partial\sigma_{\phi\rho}}{\partial\phi}+{\cot\phi\over\rho}\sigma_{\phi\rho}+{1\over\rho\sin\phi}\frac{\partial\sigma_{\theta\rho}}{\partial\theta}-{1\over\rho}\left(\sigma_{\theta\theta}+\sigma_{\phi\phi}\right)=0\\\
\theta&:\frac{\partial\sigma_{\rho\theta}}{\partial\rho}+2{\sigma_{\rho\theta}\over\rho}+{1\over\rho}\frac{\partial\sigma_{\phi\theta}}{\partial\phi}+{1\over\rho\sin\phi}\frac{\partial\sigma_{\theta\theta}}{\partial\theta}+{\sigma_{\theta\rho}\over\rho}+{\cot\phi\over\rho}\left(\sigma_{\phi\theta}+\sigma_{\theta\phi}\right)=0\\\
\phi&:\frac{\partial\sigma_{\rho\phi}}{\partial\rho}+2{\sigma_{\rho\phi}\over\rho}+{1\over\rho}\frac{\partial\sigma_{\phi\phi}}{\partial\phi}+{1\over\rho\sin\phi}\frac{\partial\sigma_{\theta\phi}}{\partial\theta}+{\sigma_{\phi\rho}\over\rho}+{\cot\phi\over\rho}\left(\sigma_{\phi\phi}-\sigma_{\theta\theta}\right)=0\end{aligned},$
(31)

and the boundary conditions corresponding to an internal pressure $p_{i}$ and
external pressure $p_{o}$ for this problem are

$\displaystyle\left.\sigma_{\rho\rho}\right|_{\rho=\rho_{i}}=-p_{i},\quad\left.\sigma_{\rho\rho}\right|_{\rho=\rho_{o}}=p_{o}.$
(32)

Specializing to the deformation (24), we conclude from the $\theta$ and $\phi$
equations that $p=p(\rho)$, and the $\rho$ equation reduces to

$\frac{\partial\sigma_{\rho\rho}}{\partial\rho}+{1\over\rho}\left(2\sigma_{\rho\rho}-\sigma_{\theta\theta}+\sigma_{\phi\phi}\right)=0.$
(33)

It follows that for any $\rho_{1}\leq\rho_{2}$,

$\left.\sigma_{\rho\rho}\right|_{\rho=\rho_{2}}-\left.\sigma_{\rho\rho}\right|_{\rho=\rho_{1}}=\int_{\rho_{1}}^{\rho_{2}}{2\over\rho}\left(\sigma_{\theta\theta}(\lambda)-\sigma_{\rho\rho}(\lambda)\right)d\rho=\int_{\rho_{1}}^{\rho_{2}}{2\over\rho}\left(\hat{\sigma}_{1}(\lambda)-\hat{\sigma}_{3}(\lambda)\right)d\rho.$
(34)

To arrive at this, we use the fact that
$\sigma_{\theta\theta}=\sigma_{\phi\phi}$ for the first equality and the
ordering of the eigenvalues of $\bm{b}$ for the second. Further, we write
$\hat{\sigma}_{i}(\lambda)=\hat{\sigma}_{i}(s(\lambda),t(\lambda))$ with a
slight abuse of notation. Note that we have to interpret this integral
appropriately by dividing the domain into various regions and using the
appropriate branch. Observe that the integrand does not involve $p$ and is
entirely constitutive. It is convenient to change variable from $\rho$ to $R$
and obtain for any $R_{1}\leq R_{2}$ that

$\left.\sigma_{\rho\rho}\right|_{R=R_{2}}-\left.\sigma_{\rho\rho}\right|_{R=R_{1}}=-\int_{R_{1}}^{R_{2}}{2\over\lambda^{3}R}\left(\hat{\sigma}_{3}(\lambda)-\hat{\sigma}_{1}(\lambda)\right)dR.$
(35)

We use this relation to analyze the three cases.

We focus on the generic Case ML assuming that $R_{i}\leq R^{*}\leq R_{o}$.
When $R^{*}\leq R\leq R_{o}$, we are in the region $L$, where the constitutive
contribution to the stress is zero. Therefore, it follows from (35) that
$\sigma_{\rho\rho}$ is constant, and therefore

$\sigma_{RR}=p_{o},\quad R^{*}\leq R\leq R_{o}.$ (36)

In particular $\sigma_{\rho\rho}(R^{*})=p_{o}$. When $R_{i}\leq R<R^{*}$, we
are in the region $M$ and we conclude from (35)

$p_{o}-p_{i}=-\int_{R_{i}}^{R^{*}}{2\over\lambda^{3}R}\left(2\hat{\sigma}_{\rho\rho}(\lambda)-\hat{\sigma}_{\theta\theta}(\lambda)\right)dR.$
(37)

We now have a system of two equations, (37) and the definition of $R^{*}$ in
(26), to solve for the two unknowns $a$ and $R^{*}$. If it turns out that
$R^{*}<R_{i}$, then we are in Case L, and if it turns out that $R^{*}>R_{i}$,
then we are in Case M.

We now specialize to two examples.

###  4.2 Inflation of a spherical balloon

We consider a spherical balloon with inner (reference) radius $R_{i}=1$ cm and
outer radius $R_{o}=1.1$ cm subjected to an internal pressure and zero
external pressure. We consider a relaxed generalized Mooney-Rivlin (RGMR)
model with parameters (16). We also consider the relaxed neo-Hookean type
Bladon-Terentjev-Warner (RBTW) model with the same $c_{1},p_{1}$. All
calculations were performed in MATLAB.

Figure 5: Normalized internal pressure vs. normalized outer radius for a
nematic balloon with anisotropy parameter $r=8$ according to the relaxed
generalized Mooney-Rivlin (RGMR) and relaxed Bladon-Terentjev-Warner (RBTW)
models. The inset in the top-right corner highlights the portion of the plot
at small $\lambda_{o}\in[1,1.3]$.

Figure 5 displays the pressure-azimuthal stretch response of the two models
when $r=8$. The entire balloon is initially in state $L$. It remains in state
$L$ for small $\lambda_{o}$ ( $\lambda_{0}\leq\lambda_{0}^{-}=1.148$, c.f.
Table 1) with zero internal pressure. As $\lambda_{o}$ increases beyond this
value, the balloon transitions to Case ML as the internal pressure increases.
The portion of the balloon in state $M$ gradually increases with increasing
stretch and pressure till $\lambda_{o}=\lambda_{0}^{+}=1.189$. The balloon is
now entirely in state $M$ (Case M), and the response stiffens. This continues
till a critical stretch at which time the response softens and eventually
suffers an instability. The behavior predicted by both models is similar up
until this point, but they diverge beyond. The RBTW model predicts complete
instability, but the RGMR model restabilizes for large $\lambda_{o}$. This is
consistent with the theoretical and experimental results on rubber balloons
[35]. Since the RBTW model is based on a Gaussian approximation to the polymer
chains, it under-predicts the response at large stretches. In summary, we
expect the nematic balloon to be initially soft as it transitions from region
$L$ to $M$, but then behave like a rubber balloon: inflating with increasing
pressure till a critical pressure, at which time it undergoes an instability
associated with a large deformation, where it restabilizes.

(a)

(b)

(c)

(d)

Figure 6: Inflation of a spherical balloon. (a) Normalize pressure vs.
normalized outer radius for varying anisotropy parameter. (b) Progression of
the various cases during inflation. (c) Expanded view of (a) for small
stretch. (d) Expanded view of (b) for small stretch.

Figures 6(a) and 6(c) show the results for varying anisotropy parameter $r$
for the RGMR material. As the pressure $p$ increases, the balloon undergoes an
instability and eventually restabilizes for all values of $r$. The response of
the balloon is stiffest in the isotropic state ($r=1$), and gets
correspondingly more compliant as $r$ increases. Further, the strain to
restabilization increases as $r$ increases. In each case, the balloon starts
in Case L, then transitions to Case ML and finally to Case M; see Figures 6(b)
and 6(d) for expanded views at small stretch values.

###  4.3 Cavitation

(a)

(b)

Figure 7: Normalized pressure vs. normalized inner radius for cavitation. (a)
Comparison of the predictions of relaxed Bladon-Terentjev-Warner (RBTW) and
relaxed generalized Mooney-Rivlin (RGMR) models for an anisotropy parameter
$r=8$. (b) Results for cavitation at varying anisotropy parameter for RGMR.

Following Gent and Lindley [18], we consider a spherical shell with an outer
radius $R_{o}=1$cm, and a very small inner radius $R_{i}=10^{-6}$ cm, and
subject it to hydrostatic tension (negative pressure or $-p_{o}$) on the
outside. We consider the same parameters as in the previous section for both
RGMR and RBTW, and the calculations were performed in MATLAB.

The hydrostatic tension (normalized by the shear modulus) is plotted as a
function of the inner radius (normalized by the reference radius) in Figure
7(a) for both the RGMR and RBTW models with an anisotropy parameter of $r=8$.
The material is in state $L$ and remains in $L$ with zero tension till
$\lambda_{i}=1.189$. It then transitions to Case ML and the pressure begins to
rise. The material begins to soften as the pressure rises, and the inner
radius becomes very large. The RBTW model predicts an asymptotic or plateau
value of pressure while the RGMR model never reaches the plateau. However,
since the normalized radius is also the azimuthal stretch and it is known that
nematic elastomers begin to tear at stretches of about 5 [38], both models
predict cavitation failure at some critical pressure.

Note that neither model predicts a transition to Case M in either case for the
stretches considered; this does not happen till $\lambda_{i}=880$; in fact
$\lambda_{i}^{+}\to\infty$ as $R_{i}/R_{o}\to 0$.

Figure 7(b) shows the cavitation results for varying anisotropy parameter $r$.
As expected, the rubber case, $r=1$, has the stiffest response, and the
response softens as $r$ increases. Therefore, we expect the nematic elastomer
to cavitate or develop internal ruptures at smaller values of the imposed
tension.

##  5 Cylindrical balloon

###  5.1 Inflation of a cylindrical shell

We study the deformation of a cylindrical shell with undeformed height $H$,
inner radius $R_{i}$, and undeformed outer radius $R_{o}$, subjected to both
internal and external pressures. We expect the deformation to be a combination
of uniform extension and radial expansion. So we make an ansatz corresponding
to Family 3 with $c=1,d=e=0,a=1/f$:

$\rho=\sqrt{{1\over f}R^{2}+b},\quad\theta=\Theta,\quad z=fZ.$ (38)

It is easy to verify that the deformation gradient and left Cauchy-Green
tensor are

$\textbf{F}=\begin{pmatrix}\frac{1}{\lambda f}&&0&&0\\\ 0&&\lambda&&\\\
0&&0&&f\end{pmatrix},\quad\textbf{b}=\begin{pmatrix}\frac{1}{\lambda^{2}f^{2}}&&0&&0\\\
0&&\lambda^{2}&&\\\
0&&0&&f^{2}\end{pmatrix},\quad\text{where}\quad\lambda={1\over R}\sqrt{{1\over
f}R^{2}+b}.$ (39)

We assume $\lambda\geq f\geq 1$ for inflation so that the principal stretches
are $\lambda$, $f$ and $1/(\lambda f)$ in the principal directions $\theta$,
$z$ and $\rho$ respectively, and $s=\lambda,t=\lambda f$.

We specialize to the case $f=1$ so that there is no axial extension. So,
$s=t=\lambda$, and the deformation at each point is a uniaxial stretch along
the azimuthal direction. We expect to see three regions (see Figure 2) with

$L=\\{R\geq R_{2}^{*}\\},\quad M=\\{R_{1}^{*}<R<R_{2}^{*}\\},\quad
S=\\{R_{1}^{*}\leq R\\}$ (40)

where

$R_{2}^{*}=b^{1/2}(r^{1/3}-1)^{-1/2},R_{1}^{*}=b/(r-1).$ (41)

Depending on the constant $b$ and the inner and outer radii, we have six
possible cases:

  * •

Case L: $R_{1}^{*}\leq R_{2}^{*}\leq R_{i}$: the entire cylinder is in region
$L$

  * •

Case ML: $R_{1}^{*}\leq R_{i}<R_{2}^{*}<R_{o}$: the inner portion of the
cylinder is in region $M$ and the outer portion is in $L$

  * •

Case M: $R_{1}^{*}\leq R_{i}<R_{o}\leq R_{2}^{*}$: the entire cylinder is in
region $M$

  * •

Case SM: $R_{i}<R_{1}^{*}<R_{o}\leq R_{2}^{*}$: the inner portion of the
cylinder is in region $S$ and the outer portion is in $M$

  * •

Case SML: $R_{i}<R_{1}^{*}<R_{2}^{*}<R_{o}$: the cylinder is in regions $S$,
then $M$, then $L$ from inside to outside

  * •

Case S: $R_{o}\leq R_{1}^{*}\leq R_{2}^{*}$: the entire cylinder is in region
$S$

A diagram illustrating the various cases is shown in Figure 8.

Figure 8: Diagram of all possible cases in the inflation of a liquid crystal
elastomer cylinder.

###  5.2 Equilibrium

The equations of equilibrium in the absence of body forces in cylindrical
coordinates are

$\begin{aligned}
\rho&:{\partial\sigma_{\rho\rho}\over\partial\rho}+{1\over\rho}{\partial\sigma_{\rho\theta}\over\partial\theta}+{\sigma_{\rho\rho}-\sigma_{\theta\theta}\over\rho}+{\partial\sigma_{\rho
z}\over\partial z}=0\\\
\theta&:{\partial\sigma_{\rho\theta}\over\partial\rho}+{1\over\rho}{\partial\sigma_{\theta\theta}\over\partial\theta}+{2\sigma_{\rho\theta}\over\rho}+{\partial\sigma_{\theta
z}\over\partial z}=0\\\ z&:{\partial\sigma_{\rho
z}\over\partial\rho}+{1\over\rho}{\partial\sigma_{\theta
z}\over\partial\theta}+{\sigma_{\rho
z}\over\rho}+{\partial\sigma_{zz}\over\partial z}=0\end{aligned}.$ (42)

The boundary conditions corresponding to an internal pressure $p_{i}$ and
external pressure $p_{o}$ for this problem are:

$\displaystyle\left.\sigma_{\rho\rho}\right|_{\rho=\rho_{i}}=-p_{i},\quad\left.\sigma_{\rho\rho}\right|_{\rho=\rho_{o}}=p_{o}.$
(43)

Specializing to the deformation (38), we conclude from the $\theta$ and $z$
equations that $p=p(\rho)$, and the $\rho$ equation reduces to

$\displaystyle\frac{d\sigma_{\rho\rho}}{d\rho}+{1\over\rho}\left(\sigma_{\rho\rho}-\sigma_{\theta\theta}\right)$
$\displaystyle=0.$ (44)

It follows that for any $\rho_{1}\leq\rho_{2}$,

$\displaystyle\left.\sigma_{\rho\rho}\right|_{\rho=\rho_{2}}-\left.\sigma_{\rho\rho}\right|_{\rho=\rho_{1}}=-\int_{\rho_{1}}^{\rho_{2}}{1\over\rho}\left(\sigma_{\rho\rho}-\sigma_{\theta\theta}\right)d\rho=-\int_{\rho_{1}}^{\rho_{2}}{1\over\rho}\left(\hat{\sigma}_{3}(\lambda)-\hat{\sigma}_{1}(\lambda)\right)d\rho,$
(45)

where we have used the ordering of the principal stretches and have written
$\hat{\sigma}_{i}(\lambda)=\hat{\sigma}_{i}(s(\lambda),t(\lambda))$ with a
slight abuse of notation. As in the case of spherical deformation, we change
variables from $\rho$ to $R$ so that

$\displaystyle\left.\sigma_{\rho\rho}\right|_{R=R_{2}}-\left.\sigma_{\rho\rho}\right|_{R=R_{1}}=-\int_{R_{1}}^{R_{2}}{1\over\lambda^{2}R}\left(\hat{\sigma}_{3}(\lambda)-\hat{\sigma}_{1}(\lambda)\right)dR.$
(46)

We use this relation in the various cases.

We focus on the generic Case SML. When $R_{2}^{*}\leq R_{0}$, we are in the
region $L$ where the constitutive contribution to the stress is zero.
Therefore, it follows from (50) that $\sigma_{\rho\rho}$ is constant, and
therefore

$\sigma_{\rho\rho}(R)=p_{o},\quad R_{2}^{*}\leq R\leq R_{o}.$ (47)

In particular, $\sigma_{\rho\rho}(R_{2}^{*})=p_{o}$. When $R_{1}^{*}\leq
R<R_{2}^{*}$, we are in the region $M$, and

$\sigma_{\rho\rho}(R)=p_{o}-\int_{R}^{R_{2}^{*}}{1\over\lambda^{2}R}\left(\hat{\sigma}^{M}_{3}(\lambda)-\hat{\sigma}^{M}_{1}(\lambda)\right)dR,\quad
R_{1}^{*}\leq R<R_{2}^{*},$ (48)

where we use the superscript $M$ to denote that we have to use the relation
for $M$ in the formula (14) for stress. In particular, we can use this formula
to compute $\sigma_{\rho\rho}(R_{1}^{*})$. Finally, when $R_{i}\leq
R<R_{1}^{*}$, we are in $S$, and

$\sigma_{\rho\rho}(R)=\sigma_{\rho\rho}(R_{1}^{*})-\int_{R}^{R_{1}^{*}}{1\over\lambda^{2}R}\left(\hat{\sigma}^{S}_{3}(\lambda)-\hat{\sigma}^{S}_{1}(\lambda)\right)dR,\quad
R_{i}\leq R<R_{1}^{*}.$ (49)

Setting $R=R_{i}$ and using the boundary condition, we obtain

$p_{o}-p_{i}=\int_{R_{i}}^{R_{1}^{*}}{1\over\lambda^{2}R}\left(\hat{\sigma}^{S}_{3}(\lambda)-\hat{\sigma}^{S}_{1}(\lambda)\right)dR+\int_{R_{1}^{*}}^{R_{2}^{*}}{1\over\lambda^{2}R}\left(\hat{\sigma}^{M}_{3}(\lambda)-\hat{\sigma}^{M}_{1}(\lambda)\right)dR.$
(50)

We solve the system of equations (50) and (40) for the unknowns
$b,R_{1}^{*},R_{2}^{*}$. Depending on the relationship between the resulting
$R_{1}^{*},R_{2}^{*}$ with $R_{i},R_{o}$, we obtain the other cases.

We illustrate this with the example of a cylindrical balloon.

###  5.3 Cylindrical balloon

We consider the inflation of a cylindrical balloon subjected to internal
pressure with no external pressure. The calculations were performed in MATLAB
with inner radius $R_{i}=1$ cm and outer radius $R_{o}=1.1$ cm, and the same
parameters as in the previous section, except $p_{2}=5$.

(a)

(b)

Figure 9: Inflation of a cylindrical balloon subject to internal pressure (a)
Normalized internal pressure vs. normalized outer radius for varying
anisotropy parameter $r$. (b) Progression of cases with increasing normalized
outer radius.

Figure 9 shows the results for the balloon inflation at varying anisotropy
parameter $r$. Let us focus on the case $r=8$, shown in magenta. The balloon
is initially in state $L$, and the pressure is zero. As internal pressure is
applied, it transitions to Case ML with an inner annulus of state $M$, and the
pressure rises with stretch. The response softens as we transition to Case M,
to a plateau before rising again (we go through Case SM and eventually Case
S). Importantly, this shows that a nematic cylinder subjected to internal
pressure can undergo a very large deformation at moderate pressures. The
response is similar but stiffer as $r$ becomes smaller.

##  6 Bending

###  6.1 Deformation

Figure 10: Schematic depicting the mid-plane of a rectangular block undergoing
bending deformation.

We study the bending of a rectangular block into an arc of a circle, a special
case of Family 1. We consider a block occupying the region
$(-W,W)\times(-L,L)\times(-H,H)$ in its reference configuration, which we
describe with rectangular Cartesian coordinates $X_{1},X_{2},X_{3}$. We
describe the deformed configuration in cylindrical coordinates $\rho,\theta,z$
and postulate the deformation

$\rho=f(X_{1}),\quad\theta=g(X_{2}),\quad z=\lambda X_{3}$ (51)

for some functions $f,g$ and constant $\lambda>1$. The covariant basis vectors
(defined as $\bm{e}_{i}=\frac{\partial\bm{x}}{\partial\xi^{i}}$, where
$\\{\xi^{i}\\}=\\{\rho,\theta,z\\}$) associated with the cylindrical
coordinate system are

$\displaystyle\bm{e}_{\rho}=\cos\theta\bm{E}_{1}+\sin\theta\bm{E}_{2},\quad\bm{e}_{\theta}=-\rho\sin\theta\bm{E}_{1}+\rho\cos\theta\bm{E}_{2}\quad\bm{e}_{z}=\bm{E}_{3},$
(52)

where $\\{{\bm{E}}_{i}\\}$ is a Cartesian frame taken to be aligned with the
reference frame for convenience. We introduce a physical (cylindrical) basis
$\bm{e}\langle i\rangle=\bm{e}_{i}/|\bm{e}_{i}|$ that is orthonormal. We can
now write the deformation gradient as

$\begin{aligned} \bm{F}=\underbrace{\frac{\partial\xi^{i}}{\partial
X^{j}}|\bm{e}_{i}||\bm{E}_{j}|^{-1}}_{\text{$F\langle
ij\rangle$}}\bm{e}\langle i\rangle\otimes\bm{E}\langle
j\rangle=\frac{df}{dX_{1}}\bm{e}\langle\rho\rangle\otimes\bm{E}\langle
1\rangle+f\frac{dg}{dX_{2}}\bm{e}\langle\theta\rangle\otimes\bm{E}\langle
2\rangle+\lambda\bm{e}\langle z\rangle\otimes\bm{E}\langle
3\rangle\end{aligned},$ (53)

and the left Cauchy-Green tensor as

$\displaystyle\bm{b}$
$\displaystyle=\underbrace{\frac{\partial\xi^{i}}{\partial
X^{j}}\frac{\partial\xi^{k}}{\partial
X^{j}}|\bm{e}_{i}||\bm{e}_{k}||\bm{E}_{j}|^{-2}}_{\text{$B\langle
ik\rangle=F\langle ij\rangle F\langle kj\rangle$}}\bm{e}\langle
i\rangle\otimes\bm{e}\langle
k\rangle=\begin{pmatrix}\left|\frac{df}{dX_{1}}\right|^{2}&0&0\\\
0&f^{2}\left|\frac{dg}{dX_{2}}\right|^{2}&0\\\ 0&0&\lambda^{2}\end{pmatrix}$
(54)

in the cylindrical physical basis.

It remains to impose incompressibility, det $\bm{F}=1$ or
$f\frac{df}{dX_{1}}\frac{dg}{dX_{2}}\lambda=1$. Now, $f$ is a function of
$X_{1}$, while $g$ is a function of $X_{2}$ alone, and $\lambda$ is a
constant. Since this identity has to hold for all $X_{1},X_{2}$, it follows
that $dg/dX_{2}$ and $fdf/DX_{1}$ have to be constant. If follows that
$g=\kappa X_{2},f^{2}=(2/(\kappa\lambda))X_{1}+\tilde{\beta}$ for constants
$\kappa$ and $\beta>0$ (we take $f\kappa>0$ without any loss of generality) so
that

$\rho=\frac{1}{\kappa}\sqrt{\frac{2\kappa}{\lambda}X_{1}+\beta},\quad\theta=\kappa
X_{2},\quad z=\lambda X_{3},$ (55)

where we set $\beta=\kappa^{2}\tilde{\beta}$. This is consistent with family 1
of Section 3 with $a=1/(\kappa\lambda),b=\kappa,c=0,d=\beta/\kappa^{2}$.

We now specialize to the case of plane strain, where $\lambda=1$. So,

$\bm{b}=\begin{pmatrix}\frac{1}{\kappa^{2}\rho^{2}}&0&0\\\
0&\kappa^{2}\rho^{2}&0\\\ 0&0&1\end{pmatrix}=\begin{pmatrix}(2\kappa
X_{1}+\beta)^{-1}&0&0\\\ 0&2\kappa X_{1}+\beta&0\\\ 0&0&1\end{pmatrix}$ (56)

in the cylindrical physical basis. The neutral axis corresponds to the surface
$X_{1}=X_{1}^{N}$, where $b_{\theta\theta}=1$ with a radius of curvature
$\rho_{N}$:

$X_{1}^{N}={1-\beta\over 2\kappa},\quad\rho_{N}={1\over\kappa};$ (57)

so $\kappa$ is the curvature. Further, we obtain the reference configuration
in the limit $\kappa\to 0,\beta\to 1$. Since $\rho>0$ and we take $\kappa>0$,
it is convenient to work in the scaled Eulerian variable

$y=\kappa\rho=\sqrt{2\kappa X_{1}+\beta}$ (58)

as the independent variable in what follows. The principal values are

$\lambda_{1}=\begin{cases}y^{-1}=\sqrt{b_{\rho\rho}}&y\leq 1,\\\
y=\sqrt{b_{\theta\theta}}&y>1\end{cases},\quad\lambda_{2}=1,\quad\lambda_{3}=\lambda_{1}^{-1}=\begin{cases}y=\sqrt{b_{\theta\theta}}&y\leq
1,\\\ y^{-1}=\sqrt{b_{\rho\rho}}&y>1\end{cases}.$ (59)

So,

$s=t=\begin{cases}y^{-1}&y\leq 1,\\\ y&y>1.\end{cases}$ (60)

Since $s=t$, the deformation is a uniaxial stretch along the azimuthal
direction. We transition from regions $L$ to $M$ when $s=t=r^{1/6}$, and from
$M$ to $S$ when $s=t=r^{1/2}$. It follows that

$\displaystyle L$ $\displaystyle=\\{r^{-1/6}\leq y\leq r^{1/6}\\},$ (61)
$\displaystyle M$ $\displaystyle=\\{r^{-1/2}\leq
y<r^{-1/6}\\}\cup\\{r^{1/6}<X_{1}\leq r^{1/2}\\},$ (62) $\displaystyle S$
$\displaystyle=\\{y<r^{-1/2}\\}\cup\\{r^{-1/2}<y\\}.$ (63)

Defining

$y^{\pm}=\sqrt{\pm 2\kappa
W+\beta},\quad\mbox{and}\quad\rho^{\pm}=\frac{1}{\kappa}\sqrt{+2\kappa
W+\beta},$ (64)

we have nine cases666with the terminology SMLMS denoting that beam goes
through regions $S,M,L,M,S$ from bottom ($X_{1}=-W$) to top ($X_{1}=W$), etc.:

  * •

Case L: $r^{-1/6}\leq y_{W}^{-}<y_{W}^{+}\leq r^{1/6}$, i.e.,
$X_{LM}^{-}\leq-W<W\leq X_{LM}^{+}$;

  * •

Case ML: $r^{-1/2}\leq y_{W}^{-}<r^{-1/6}<y_{W}^{+}\leq r^{1/6}$, i.e.,
$X_{MS}^{-}\leq-W<X_{LM}^{-}<W\leq X_{LM}^{+}$;

  * •

Case LM: $r^{-1/6}\leq y_{W}^{-}<r^{1/6}<y_{W}^{+}\leq r^{1/2}$, i.e.,
$X_{LM}^{-}<-W<X_{LM}^{+}\leq W\leq X_{MS}^{+}$;

  * •

Case SML: $y_{W}^{-}\leq r^{-1/2}<y_{W}^{+}\leq r^{1/6}$, i.e., $-W\leq
X_{SM}^{-}<W<X_{LM}^{+}$;

  * •

Case LMS: $r^{-1/6}\leq y_{W}^{-}\leq r^{1/6}\leq r^{1/2}<y_{W}^{+}$, i.e.,
$X_{LM}^{-}\leq-W\leq X_{LM}^{+}\leq X_{MS}^{+}<W$;

  * •

Case MLM: $r^{-1/2}\leq y_{W}^{-}<r^{-1/6}<r^{1/6}<y_{W}^{+}\leq r^{-1/2}$,
i.e., $X_{MS}^{-}\leq-W<X_{LM}^{-}<X_{LM}^{+}<W\leq X_{MS}^{+}$;

  * •

Case SMLM: $y_{W}^{-}<r^{-1/2}<r^{1/6}<y_{W}^{+}\leq r^{1/2}$, i.e.,
$-W<X_{MS}^{-}<X_{LM}^{+}<W\leq X_{MS}^{+}$;

  * •

Case MLMS: $r^{-1/2}\leq y_{W}^{-}<r^{-1/6}<r^{1/2}<y_{W}^{+}$, i.e.,
$X_{MS}^{-}\leq-W<X_{LM}^{-}<X^{+}_{MS}<W$;

  * •

Case SMLMS: $y_{W}^{-}<r^{-1/2}<r^{1/2}<y_{W}^{+}$, i.e.,
$-W<X_{MS}^{-}<X_{MS}^{+}<W$,

where

$X_{LM}^{\pm}={r^{\pm 1/3}-\beta\over 2\kappa},\quad X_{MS}^{\pm}={r^{\pm
1}-\beta\over 2\kappa}.$ (65)

The various regions are highlighted in the $\kappa-\beta$ plane in Figure 11.
We only label five regions since we shall see that the other regions do not
satisfy the equilibrium equations.

Figure 11: The various cases in the $\kappa-\beta$ plane.

###  6.2 Equilibrium

We work in the current configuration in cylindrical coordinates. The equations
of equilibrium in the absence of body forces are (42), as in the case of the
inflation of the cylinder. We are given that the faces $X_{1}=\pm W$ are
traction-free so that

$\displaystyle\sigma_{\rho\rho}\big{\rvert}_{\rho=\rho^{\pm}}=0,\quad\rho^{\pm}=\kappa
y^{\pm}_{W}.$ (66)

The total axial force per unit depth $F$ is zero, while the total bending
moment per unit depth is given as $M$:

$F=\int_{\rho^{-}}^{\rho^{+}}\sigma_{\theta\theta}d\rho=0,\quad
M=\int_{\rho^{-}}^{\rho^{+}}\rho\sigma_{\theta\theta}d\rho.$ (67)

Note that the equilibrium of the entire domain and (66) automatically implies
$F=0$.

As in the previous sections, the $\theta$ and $z$ components of the
equilibrium equations lead us to conclude that the Lagrange multiplier
$p=p(\rho)$, and the $\rho$ component reduces to (44). Recalling that
$\sigma_{\rho\rho}=-p+\hat{\sigma}_{\rho\rho}$ and
$\sigma_{\theta\theta}=-p+\hat{\sigma}_{\theta\theta}$, we can rewrite the
$\rho$ equation (44) as

${dp\over
d\rho}=\frac{d\hat{\sigma}_{\rho\rho}}{d\rho}+{1\over\rho}\left(\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta}\right).$
(68)

Note that all the terms on the right-hand side are constitutively determined,
and therefore known up to the constants $\kappa,\beta$. We can integrate this
to obtain the pressure distribution

$p(\rho)=p_{0}+\int_{\rho^{-}}^{\rho}\left(\frac{d\hat{\sigma}_{\rho\rho}}{d\rho}+{1\over
s}\left(\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta}\right)\right)ds,$
(69)

where $p_{0}=p(\rho^{-})$. Note again that care should be taken to divide the
domain into the different regions and use the appropriate branch of the
constitutive relation in these regions. We may use (66) to infer that

$p_{0}=\hat{\sigma}_{\rho\rho}(\rho^{-}).$ (70)

It follows that

$p(\rho)=\hat{\sigma}_{\rho\rho}(\rho^{-})+\int_{\rho^{-}}^{\rho}\left(\frac{d\hat{\sigma}_{\rho\rho}}{d\rho}+{1\over
s}\left(\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta}\right)\right)ds.$
(71)

Finally, we can integrate the equilibrium equation (44) from $\rho=\rho^{-}$
to $\rho=\rho^{+}$ and apply the boundary condition (66) to obtain

$0=\int_{\rho^{-}}^{\rho^{+}}\frac{1}{\rho}(\sigma_{\rho\rho}-\sigma_{\theta\theta})d\rho=\int_{\rho^{-}}^{\rho^{+}}\frac{1}{\rho}(\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta})d\rho=\int_{y^{-}_{W}}^{y^{+}_{W}}\frac{1}{y}(\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta})dy.$
(72)

Given $\kappa$, we solve (72) for $\beta$; we can now determine the moment $M$
from (67)2:

$\displaystyle M$ $\displaystyle=$
$\displaystyle\int_{\rho^{-}}^{\rho^{+}}\rho\sigma_{\theta\theta}d\rho=\int_{\rho^{-}}^{\rho^{+}}\rho(-p(\rho)+\hat{\sigma}_{\theta\theta}(\rho))d\rho$
(73) $\displaystyle=$
$\displaystyle\int_{\rho^{-}}^{\rho^{+}}\rho\left(-\hat{\sigma}_{\rho\rho}(\rho^{-})-\int_{\rho^{-}}^{\rho}\left(\frac{d\hat{\sigma}_{\rho\rho}}{d\rho}+{1\over
s}\left(\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta}\right)\right)ds+\hat{\sigma}_{\theta\theta}(\rho)\right)d\rho.$
(74)

Note

$\int_{\rho^{-}}^{\rho^{+}}\rho\left(\int_{\rho^{-}}^{\rho}\frac{d\hat{\sigma}_{\rho\rho}}{d\rho}(s)ds\right)d\rho=\int_{\rho^{-}}^{\rho^{+}}\rho(\hat{\sigma}_{\rho\rho}(\rho)-\hat{\sigma}_{\rho\rho}(\rho^{-}))d\rho,$
(75)

and

$\int_{\rho^{-}}^{\rho^{+}}\rho\left(\int_{\rho^{-}}^{\rho}{1\over
s}\left(\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta}\right)ds\right)d\rho=\frac{(\rho^{+})^{2}}{2}\int_{\rho^{-}}^{\rho^{+}}\frac{1}{\rho}(\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta})d\rho-\int_{\rho^{-}}^{\rho^{+}}\frac{\rho}{2}(\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta})d\rho$
(76)

using integration by parts. The first term on the right-hand side above is
zero by (72). So,

$M=\frac{1}{2}\int_{\rho^{-}}^{\rho^{+}}\rho(\hat{\sigma}_{\theta\theta}-\hat{\sigma}_{\rho\rho})d\rho=\frac{1}{2\kappa^{2}}\int_{y^{-}_{W}}^{y^{+}_{W}}y(\hat{\sigma}_{\theta\theta}-\hat{\sigma}_{\rho\rho})dy.$
(77)

Note that in evaluating both (72) and (77), we need to divide the domain into
the various regions and use the appropriate branch of the constitutive
relation.

Now, note that in light of (59) and our constitutive relation,

$\hat{\sigma}_{\rho\rho}-\hat{\sigma}_{\theta\theta}\begin{cases}<0\quad&y<r^{-1/6}\\\
=0&r^{-1/6}\leq y\leq r^{1/6}\\\ >0&y>r^{1/6}\end{cases}.$ (78)

It follows therefore that the cases ML, LM, SML and LMS cannot satisfy the
equilibrium condition (72). It also follows that we are in the Case L for
small $\kappa$ and remain there till we simultaneously satisfy the conditions
$y^{-}_{W}=r^{-1/6}$ and $y^{+}_{W}=r^{1/6}$. In other words,

$M=0,\quad
0\leq\kappa\leq\kappa_{0}:=\frac{1}{4W}\left(r^{1/3}-r^{-1/3}\right).$ (79)

At this point, $\beta=(r^{1/3}+r^{-1/3})/2$. As $\kappa$ increases beyond
$\kappa_{0}$, we have Cases MLM, SMLM/MLMS and possibly SMLMS, and the moment
$M$ is positive.

Figure 12: Progression of the bending solution through individual case
numbers.

###  6.3 Case study

The stresses, forces, and moments as well as the deformed radii were all
solved for using Mathematica using the following parameters: $c_{1}=1.03\cdot
10^{5}$ Pa, $c_{2}=1.96\cdot 10^{2}$ Pa, $d_{1}=1.63\cdot 10^{-2}$ Pa,
$p_{1}=1$, $p_{2}=5$, and $q_{1}=1$. The results are shown in Figure 12. We
find in this example that we observe only three cases: L, MLM and SMLMS.

##  7 Acknowledgments

We are grateful for the financial support of the US Air Force Office of
Scientific Research through the MURI Grant No. FA9550-16-1-0566.

## References

  * [1] V. Agostiniani and A. DeSimone.  Ogden-type energies for nematic elastomers.  International Journal of Non-Linear Mechanics, 47:402–412, 2012\.
  * [2] H. Alexander.  Tensile instability of initially spherical balloons.  International Journal of Engineering Science, 9:151–160, 1971.
  * [3] C. P. Ambulo, J. J. Burroughs, J. M. Boothby, H. Kim, M. R. Shankar, and T. H. Ware.  Four-dimensional Printing of Liquid Crystal Elastomers.  ACS Applied Material Interfaces, 9:37332–37339, 2017.
  * [4] G. Baardink and P. Cesana.  Torsion of a liquid crystal elastomer cylinder.  Working notes: private communication, 2019.
  * [5] J. M. Ball.  Discontinuous equilibrium solutions and cavitation in nonlinear elasticity.  Philosophical Transactions of the Royal Society of London A, 306:557–611, 1982.
  * [6] J. S. Biggins, E. M. Terentjev, and M. Warner.  Semisoft response of nematic elastomers to complex deformations.  Physical Review E, 78:041704, 2008.
  * [7] J. S. Biggins, M. Warner, and K. Bhattacharya.  Supersoft elasticity in polydomain nematic elastomers.  Physical Review Letters, 103, 2009.
  * [8] P. Bladon, M. Terentjev, and M. Warner.  Transitions and instabilities in liquid-crystal elastomers.  Physical Review E, 47:3838–3840, 1993.
  * [9] P. Cesana, P. Plucinsky, and K. Bhattacharya.  Effective Behavior of Nematic Elastomer Membranes.  Archive for Rational Mechanics and Analysis, 218:863–905, 2015\.
  * [10] S. Conti, A. Desimone, and G. Dolzmann.  Semisoft elasticity and director reorientation in stretched sheets of nematic elastomers.  Physical Review E, 66:061710, 2002.
  * [11] S. Conti, A. DeSimone, and G. Dolzmann.  Soft elastic response of stretched sheets of nematic elastomers: a numerical study.  Journal of the Mechanics and Physics of Solids, 50:1431–1451, 2002\.
  * [12] B. Dacorogna.  Direct methods in the calculus of variations.  Springer, New York, 2008.
  * [13] P.-G. de Gennes.  Réflexions sur un type de polymères nématiques.  Comptes rendus de l’Académie des Sciences, Série B, 281:101–103, 1975.
  * [14] A. DeSimone and G. Dolzmann.  Material instabilities in nematic elastomers.  Physica D, 136:175–191, 2000.
  * [15] A. DeSimone and G. Dolzmann.  Macroscopic Response of Nematic Elastomers via Relaxation of a Class of SO(3)-Invariant Energies.  Archive for Rational Mechanics and Analysis, 161:181–204, 2002\.
  * [16] J. L. Ericksen.  Deformations possible in every isotropic, incompressible, perfectly elastic body.  Zeitschrift für angewandte Mathematik und Physik ZAMP, 5:466–489, 1954.
  * [17] H. Finkelmann, H.-J. Kock, and G. Rehage.  Investigations on liquid crystalline polysiloxanes. 3. liquid crystalline elastomers – a new type of liquid crystalline material.  Macromolecules. Rapid Communications, 2:317–322, 1981.
  * [18] A. N. Gent and P. B. Lindley.  Internal rupture of bonded rubber cylinders in tension.  Proceedings of the Royal Society of London A, 249:195–205, 1959\.
  * [19] A. Giudici and J. S. Biggins.  Giant deformations and soft-inflation in LCE balloons.  Europhysics Letters, 132:36001, 2020.
  * [20] Q. He, Y. Zheng, Z. Wang, X. He, and S. Cai.  Anomalous inflation of a nematic balloon.  Journal of the Mechanics and Physics of Solids, 142:104013, Sept. 2020.
  * [21] W. W. Klingbeil and R. T. Shield.  On a class of solutions in plane finite elasticity.  Zeitschrift für angewandte Mathematik und Physik ZAMP, 17:489–511, 1966.
  * [22] I. Kundler and H. Finkelmam.  Strain-induced director reorientation in nematic liquid single crystal elastomers.  Macromol. Rapid Commun, 16:679–686, 1995.
  * [23] J. Küpfer and H. Finkelmann.  Nematic liquid single crystal elastomers.  Macromolecular Rapid Communications, 12:717–726, 1991.
  * [24] L. Lanzoni and A. M. Tarantino.  The Bending of Beams in Finite Elasticity.  Journal of Elasticity, 139:91–121, 2020.
  * [25] V. Lee and K. Bhattacharya.  Actuation of cylindrical nematic elastomer balloons.  Journal of Applied Physics, 129:114701, 2021.
  * [26] L. A. Mihai and A. Goriely.  Instabilities in liquid crystal elastomers.  MRS Bulletin, 46:784–794, 2021.
  * [27] L. Petridis and E. Terentjev.  Nematic-isotropic transition with quenched disorder.  Physical Review E, 74:051707, 2006.
  * [28] R. S. Rivlin.  Large elastic deformations of isotropic materials VI. Further results in the theory of torsion, shear and flexure.  Philosophical Transactions of the Royal Society of London A, 845:173–195, 1949.
  * [29] R. S. Rivlin and E. K. Rideal.  Large elastic deformations of isotropic materials. III. Some simple problems in cyclindrical polar co-ordinates.  Philosophical Transactions of the Royal Society of London A, 240:509–525, 1948.
  * [30] R. S. Rivlin and E. K. Rideal.  Large elastic deformations of isotropic materials. V. The problem of flexure.  Proceedings of the Royal Society of London A, 195:463–473, 1949\.
  * [31] R. S. Rivlin and E. K. Rideal.  Large elastic deformations of isotropic materials VI. Further results in the theory of torsion, shear and flexure.  Philosophical Transactions of the Royal Society of London A, 242:173–195, 1949.
  * [32] M. O. Saed, W. Elmadih, A. Terentjev, D. Chronopoulos, D. Williamson, and E. M. Terentjev.  Impact damping and vibration attenuation in nematic liquid crystal elastomers.  Nature Communications, 12:6676, 2021.
  * [33] M. Singh and A. C. Pipkin.  Note on Ericksen’s problem.  Zeitschrift für angewandte Mathematik und Physik ZAMP, 16(5):706–709, 1965.
  * [34] H. Tokumoto, H. Zhou, A. Takebe, K. Kamitani, K. Kojio, A. Takahara, K. Bhattacharya, and K. Urayama.  Probing the in-plane liquid-like behavior of liquid crystal elastomers.  Science Advances, 7:abe9495, 2021.
  * [35] L. R. G. Treloar.  The Physics of Rubber Elasticity.  Oxford University Press, Oxford, 3rd edition, 2005.
  * [36] K. Urayama, E. Kohmon, M. Kojima, and T. Takigawa.  Polydomain\Monodomain Transition of Randomly Disordered Nematic Elastomers with Different Cross-Linking Histories.  Macromolecules, 42:4084–4089, May 2009.
  * [37] T. H. Ware, M. E. McConney, J. J. Wie, V. P. Tondiglia, and T. J. White.  Voxelated liquid crystal elastomers.  Science, 347:982–984, 2015.
  * [38] M. Warner and E. M. Terentjev.  Liquid Crystal Elastomers.  Oxford University Press, Oxford, 2003.
  * [39] T. J. White and D. J. Broer.  Programmable and adaptive mechanics with liquid crystal polymer networks and elastomers.  Nature Materials, 14:1087–1098, 2015.
  * [40] C. M. Yakacki, M. Saed, D. P. Nair, T. Gong, S. M. Reed, and C. N. Bowman.  Tailorable and programmable liquid-crystalline elastomers using a two-stage thiol-acrylate reaction.  RSC Advances, 5:18997–19001, 2015.
  * [41] A. Yavari.  Universal deformations in inhomogeneous isotropic nonlinear elastic solids.  Proceedings of the Royal Society A: Mathematical, Physical and Engineering Sciences, 477:20210547, 2021.
  * [42] H. Zhou and K. Bhattacharya.  Accelerated computational micromechanics and its applications to nematic elastomers.  Journal of the Mechanics and Physics of Solids, 153:104470, 2021\.
  * [43] L. Zubov and M. Karyakin.  Nonlinear deformations of a cylindrical pipe with pre-stressed thin coatings.  Mathematics and Mechanics of Solids, 27:1703–1720, 2022.