Nematic elastomers are programmable soft materials that display large, reversible, and predictable deformation under an external stimulus such as a change in temperature or light. While much of the work in the field has focused on actuation from flat sheets, recent advances in 3D printing and other methods of directed synthesis have motivated the study of actuation of curved shells. Snap-through buckling has been a topic of particular interest. In this work, we present theoretical calculations to motivate another mode of actuation that combines programmable soft materials as well as instabilities associated with large deformation. Specifically, we analyze the deformation of a cylindrical shell of a patterned nematic elastomer under pressure, show that it can undergo an enormous change of volume with changing temperature and suggest its application as a pump with extremely high ejection fraction.

## I. INTRODUCTION

Liquid crystal elastomers (LCEs) are remarkable stimuli-responsive materials that have recently been explored for their use in actuation.^{29} LCEs are elastomers with liquid crystal mesogens incorporated into the underlying polymer chains. Of particular interest are nematic elastomers that undergo an isotropic-to-nematic transition accompanied by a significant stretch (by as much as a factor of two or more) along the nematic director and lateral contraction, as shown schematically in Fig. 1. The stiff rod-like nematic mesogens are randomly oriented at high temperature (isotropic state) but align along a particular direction due to steric interactions below a critical temperature (the nematic state). This phase transition induces a deformation of the underlying polymer chains in the nematic elastomer. The degree of orientation, and consequently the change of shape, depend on temperature.

Nematic elastomers have been exploited for programmable actuation and shape-morphing of thin sheets. Modes *et al.*^{17,18} suggested that if sheets of nematic elastomer with prescribed director patterns were fabricated in the nematic state and subsequently heated, they could deform out of plane into three-dimensional shapes. For example, a +1 disclination with an azimuthal director distribution would deform into a cone. This was demonstrated in nematic glasses by de Haan *et al.*^{7} Ware *et al.*^{28} developed a method of synthesizing nematic elastomers where the director pattern could be written pixel by pixel on flat sheets and demonstrated the formation of these cones. Moreover, they showed that this actuation was extremely robust, as the cone-lifting weights were many hundreds of times larger than the structure itself. Since then, there have been a number of other studies on nematic elastomers,^{2,16,19,22,23} including the inverse problem of identifying the director pattern that would lead to a given actuated shape.^{3,23} All of these works address the programming and actuation of initially flat sheets.

Recent advances in 3D printing and other methods of directed synthesis have enabled the synthesis of curved shells,^{4,12} and such structures change shape upon heating. In particular, Ambulo *et al.* demonstrated dramatic snap-through buckling of structures involving regions of positive Gauss curvature and regions of negative Gauss curvature.^{4} More recently, magnetic fields have been used to independently control director orientation during 3D printing.^{26} These developments in synthesis techniques motivate the current work.

The goal of this work is twofold. The first is to explore the combination of programmed synthesis of nematic shells and the geometric instabilities associated with the large deformation of slender structures. Similar instabilities have been exploited in other stimuli-responsive materials including electroactive materials.^{13,25,32} In this work, we focus on the so-called aneurysm instability of pressurized cylinders.^{10} As observed in long toy balloons, one observes a discontinuous change of radius (or volume) with an increase of pressure: typically, the balloon inflates until it reaches a particular radius, beyond which point a bump (aneurysm) with a significantly larger radius appears in this region, and it propagates through the entire balloon before the radius further increases. We explore the response of a cylindrical shell made of a nematic elastomer and study how the isotropic–nematic phase transition affects this instability. Our work is closely related to that of Giudici and Biggins,^{11} who recently studied the ballooning instability in both nematic and isotropic LCEs using a Gent-style energy. He *et al.*^{14} have studied the anomalous behavior of (isotropic-genesis polydomain) nematic balloons under tension. We then show how this instability can be used as a high ejection-fraction pump. The second goal is to study actuation and shape-morphing in the presence of mechanical loads. The prior literature has largely focused on free recovery.

## II. LARGE DEFORMATION MODEL OF NEMATIC ELASTOMERS

We begin with the neo-classical theory of nematic elastomers following Bladon *et al*.^{6,29} The state of a liquid crystal elastomer is described by an anisotropy parameter $r$, a director $n$, and the deformation gradient $F$ relative to a stress-free reference configuration with anisotropy parameter $r0$ and director $n0$. The anisotropy parameter is a function of temperature with $r=1$ in the isotropic state above the transformation and gradually increases with decreasing temperature so that $r>1$ in the nematic state. We consider the material to be incompressible so that $detF=1$. The neo-classical theory considers the entropy of the polymer chains in the Gaussian approximation, and the free energy density is given as

where $\mu $ is the shear modulus of the material and

are the step-length tensors in the current and reference configurations that collect the anisotropy parameter and the director. It is easy to show that

where $WNH(F)=\mu 2(trC(F))\u22123)$, with $C(F)=FTF$ being the neo-Hookean energy density, which describes the entropy of polymer chains in ordinary rubber in the Gaussian approximation.^{27}

The neo-classical theory is known to describe complex features of nematic elastomers at finite, but moderate, deformation. However, at extremely large stretches, the Gaussian approximation does not hold, and this theory does not adequately describe the stiffening much like its neo-Hookean counterpart. Various constitutive relations are used to describe rubber in this high-stretch regime. A common feature of many of these models is that the energy density depends only on principle stretches $\lambda i$ of $F$ [equivalently the eigenvalues $\lambda i2$ of $C(F)$],

For example, in the Ogden model,^{20} the energy density is

where $N$, $\mu p$, and $\beta p$ are material constants. The shear modulus is $\mu =12\u2211p=1N\mu p\beta p$. When $N=1$ and $\beta 1=2$, the Ogden energy is the neo-Hookean energy, and when $N=2$, $\beta 1=2$, and $\beta 2=\u22122$, the Ogden energy is the Mooney–Rivlin energy. We use the Ogden energy to demonstrate our results following Ref. 25, though we can adapt them to any constitutive relation that describes the high-stretch behavior. We adopt the elastic energy density (5) to nematic elastomers analogously to (4). See Ref. 1 for similar energies and their relaxation in the ideal case. Other approaches have been proposed to capture the high-stretch regime including the logarithmic correction by Gent,^{9} which was used by Giudici and Biggins^{11} in their work.

Furthermore, the cross-link density and the polymer network may carry an imprint of the initial director, and this leads to breaking of symmetry (isotropy), leading to the preference of the director to remain in the original orientation. Such an interaction can be described using an additional non-ideal energy density,^{5}

Note that this energy is minimized when $n=n0$. Putting these together, we take the energy density of the nematic elastomer to be

For future use, we note a particular invariance of this energy density. Let $Q$ be a rotation tensor that leaves the reference director invariant, $Qn0=\xb1n0$. Then, we claim that

Note that

where we have used the invariance of $n0$ under $Q$ in the third equality. It follows that both tensors have the same eigenvalues and thus the same Ogden energy density. A similar calculation holds for the non-ideal energy density as well, thereby establishing (9).

## III. INFLATION OF A NEMATIC CYLINDER

Consider a cylindrical shell of initial (reference) length $H$, inner radius $Ri$, and undeformed outer radius $Ro$ subjected to an internal pressure $p$. Following Rivlin^{24} and Ericksen,^{8} the deformation of the cylinder is described by the universal volume-preserving deformation involving radial expansion, axial extension, and torsion (see Fig. 2). The mapping is

where ${R,\Theta ,Z}$ and ${\rho ,\theta ,z}$ denote the cylindrical coordinate system in the reference and deformed coordinate systems, respectively. $\rho (R)$ describes the radial expansion, $D$ is the twist, and $\xi $ is the axial stretch. The deformation gradient in the cylindrical coordinate system is

where we have introduced the hoop stretch $\lambda (R)=\rho (R)/R$ and used the incompressibility to obtain the second equality. There is an off-diagonal term in the deformation gradient because we would like to allow shear or twist that may accompany director reorientation. We will see later that this indeed plays a role.

We assume that directors both in the reference and deformed configuration are tangential to the cylinder and make an angle $\varphi 0$ and $\varphi $, respectively, with the azimuthal coordinate. Thus, in cylindrical coordinates,

The total potential energy of the system is

where $\Delta V$ is the difference in the deformed and undeformed volumes. Applied to a balloon with height $H$, we obtain

Above we have assumed that the shell is thin, $T:=(Ro\u2212Ri)\u226aRi$, to evaluate the integral.

For a given pressure $p$ and anisotropy parameter $r$, we can now find the equilibrium as

Physically, these equations describe the balance between the hoop stress in the cylinder and the internal pressure, the balance between the axial stress and the internal pressure, the balance of torque, and the balance of internal (material) torque on the director, respectively.

To demonstrate the results, we consider a cylinder where the initial director is axial ($\varphi 0=90\xb0$) and which is mildly nematic with initial anisotropy parameter $r0=2$. The rest of the parameters are shown in Table I and were chosen to be broadly consistent with an experiment conducted in our laboratory. We fix the current anisotropy parameter $r$ and the hoop stretch $\lambda $ and solve (16) for the pressure $p$, the axial stretch $\xi $, the twist $D$, and the current director angle $\varphi $. We find that the system has two solutions, shown in Fig. 3 for four different current anisotropy parameters *r* and the classical neo-Hookean case *r*_{0} = *r* = 1. There is a third unstable solution where the director does not rotate, which we ignore. Note that the pressure has been normalized as $p^=pRo\mu (Ro\u2212Ri)$, where $\mu =12\u2211p=13\mu p\beta p=6.80\xd7104$ Pa represents the reference shear modulus.

Inner radius | R_{i} | 1 cm |

Outer radius | R_{o} | 1.05 cm |

Height of cylinder | H | 5 cm |

Initial director angle | ϕ_{0} | 90° |

Initial anisotropy parameter | r_{0} | 2 |

Non-ideality parameter | α | 0.3 |

Ogden model shear modulus | μ_{1} | 1.0 × 10^{5} Pa |

Ogden model shear modulus | μ_{2} | 1.904 762 × 10^{2} Pa |

Ogden model shear modulus | μ_{3} | −1.587 3 × 10^{3} Pa |

Ogden model constant | β_{1} | 1.3 |

Ogden model constant | β_{2} | 6 |

Ogden model constant | β_{3} | −3 |

Inner radius | R_{i} | 1 cm |

Outer radius | R_{o} | 1.05 cm |

Height of cylinder | H | 5 cm |

Initial director angle | ϕ_{0} | 90° |

Initial anisotropy parameter | r_{0} | 2 |

Non-ideality parameter | α | 0.3 |

Ogden model shear modulus | μ_{1} | 1.0 × 10^{5} Pa |

Ogden model shear modulus | μ_{2} | 1.904 762 × 10^{2} Pa |

Ogden model shear modulus | μ_{3} | −1.587 3 × 10^{3} Pa |

Ogden model constant | β_{1} | 1.3 |

Ogden model constant | β_{2} | 6 |

Ogden model constant | β_{3} | −3 |

Consider the first solution, Figs. 3(a)–3(d). We observe that for $r0\u2260r$, hoop stretch vs pressure does not pass through $(1,0)$ but through $((r/r0)\u22121/6,0)$ since the change of the anisotropy parameter gives rise to a spontaneous deformation of the cylinder. The hoop stretch vs pressure is non-monotone [Fig. 3(a)]: the pressure initially increases but then drops before increasing again with increasing hoop stretch. This reflects the well-known balloon instability: with increasing pressure, the radius increases until it reaches a critical pressure at which it jumps to a large radius. The onset and the extent of this instability are amplified in nematic elastomers due to the rotation of the director. Figure 3(c) shows that the director begins to rotate with inflation, reaching the hoop direction asymptotically. To understand this, an increase in radius increases the volume more than an increase in the axial stretch, since the former leads to an increase of included area rather than length. Therefore, the pressure seeks to increase the circumference by reorienting the director. This reorientation also leads to a decrease of the axial stretch [Fig. 3(b)]. Consequently, the axial stretch is also non-monotone: it decreases during reorientation but increases again as the director stabilizes. Finally, the reorientation leads to a twist in the cylinder [Fig. 3(d)]. The magnitude of all of these trends increase with increasing anisotropy parameter $r$. In particular, the critical pressure decreases and the change of hoop stretch increases with increasing anisotropy parameter $r$.

The reorientation, however, is resisted by the non-ideality, as shown by varying the non-ideality parameter $\alpha $ in Fig. 4. Note that the director rotation from the vicinity of the initial orientation $\varphi 0=\pi 2$ is increasingly delayed as the non-ideality parameter increases. The balance between the pressure-assisted reorientation and the non-ideality-mediated resistance leads to the observed behavior.

The second solution, Figs. 3(a), 3(b), 3(e), and 3(f), is very similar to the first, except that the reorientation and twist change sign. The pressure vs hoop-stretch and the axial stretch vs hoop-stretch curves remain unchanged. Consequently, both solutions have the same pressure vs volume curves, which are shown in Fig. 5(a). The volume strain is plotted on a logarithmic scale due to the dramatic change of volume during the instability.

It is useful to understand the origin of the two solutions. The material is not chiral, and neither is the initial configuration. Therefore, the breaking of the chiral symmetry by rotation of the director has to be accompanied by a symmetry-related counterpart. To elaborate on this, recall the invariance (9). Let $Q$ be a $180\xb0$ rotation about the azimuthal direction,

Note that $Qn0=\u2212n0$ so that it satisfies the requirement for (9). It is easy to check that for $F$ and $n$ in (12) and (13),

Thus, the invariance (9) implies that any solution with chirality has a symmetric counterpart with the same radial and azimuthal stretches.

The presence of the two symmetric solutions enables the formation of stripe domains that avoid overall torsion, as shown in Fig. 5(b). We divide the cylinder into short cylindrical rings and alternate between the two solutions. This leads to a continuous deformation, where one ring twists one way and the other ring the other way in an alternating pattern, but they meet continuously across the boundaries as indicated by the initially straight fiducial dashed line shown in the figure. The overall torsion is zero, while the overall hoop and axial stretch are as before, leading to the pressure–volume curve shown in Fig. 5(a), where volume strain is defined as the current volume divided by the reference volume.

Stripe domains are widely observed in nematic elastomers, especially in uniaxial tension, where rigid grips prevent any shear.^{29} In uniaxial tension of a nematic sheet along an axis that is perpendicular to the initial director orientation, director rotation accommodates stretch but causes shear. However, shear breaks the symmetry and, therefore, there are two solutions (rotation to the right or left), which alternate to form the stripe domains. The domains are fine, typically with the width of micrometers, and the interfaces are very sharp, with a width of nanometers. The stripe domains in Fig. 5(b) are the exact analogs of those in uniaxial tension.

## IV. PUMP

The pressure–volume curves in Fig. 5(a) motivate the application of this cylindrical nematic elastomer balloon as a pump. Recall that the anisotropy parameter $r$ depends on temperature, and, therefore, the four pressure–volume curves represent four distinct temperatures. In a typical monodomain nematic elastomer, $r=2$ at a high temperature of about $85\xb0C$, while $r=8$ at a low temperature of about $25\xb0C$.^{29} These two pressure–volume curves are re-plotted in Fig. 6 as the hot and cold nematic elastomers. An important observation is that the lower-critical pressure (point E) of the pressure–volume curve at the high temperature is higher than the upper-critical pressure (point B) of the pressure–volume curve at the low temperature. This enables the operation as a pump between an inlet pressure $pi$ and outlet pressure $po$, where $pB<pi\u2264po\u2264pE$.

Also shown in the figure are the isotherms (pressure–volume relation) of a fixed mass of fluid, in this case air, at the hot and the cold temperatures of $85\xb0C$ and $25\xb0C$, respectively. These isotherms were calculated using the ideal gas law. The cold fluid isotherm is given by $pcold(V)=nCRTcold/V$, where $nC$ is the number of moles of air at point C, $R=8.3145$ is the ideal gas law constant, and $Tcold$ is the cold fluid temperature. The hot fluid isotherm pressure is analogously given by the corresponding number of moles at point D and the hot fluid temperature: $phot(V)=nDRThot/V$.

The pump operates as follows. Let us begin at the high temperature with the outlet closed and the inlet open so that the nematic pump is at the point A. Now, cool the pump with the inlet open so that the pressure remains at $pi$. On cooling fully to the cold temperature when $r=8$, the only equilibrium solution is point C, which has a very large volume. So, the balloon would draw in a large volume of air from A to C. We note that this process does not proceed smoothly. As the pump is cooled from the high temperature, the volume changes gradually until the temperature when the upper-critical pressure (the point corresponding to B at the intermediate temperature) equals $pi$. At this point, there will be an instability (likely accompanied by an aneurysm), and the volume jumps to something close to C. This instability has been analyzed by Giudici and Biggins^{11} and may be of interest in microfluidics. Subsequent cooling takes it to point C.

Now, close the inlet and start heating the pump. The mass of fluid in the pump is fixed, and so its behavior shifts from that of the cold isotherm to that of the hot isotherm. In the interim, the pressure–volume curve of the pump also changes to that of the hot material. Therefore, the equilibrium shifts from C to D. Now, open the outlet so that the pressure decreases to $po$. The only available state in the hot pump is at F, and so the pump goes from point D, with very large volume, to point F, with small volume, expelling the fluid. This is again accompanied by an instability from E to F. Closing the outlet and opening the inlet takes us from F to A, resetting the pump.

A pump can be characterized by its ejection fraction. In this case, the ejection fraction is

which means that $98.6%$ of the fluid is pumped out of a filled balloon during each cycle. This is extremely high: a normal human heart has a left ventricular ejection fraction between $50%$ and $70%$, according to the American College of Cardiology.^{15} A plot of ejection fraction for fixed $r0=2$ and varying $r$ can be seen in Fig. 7.

## V. CONCLUSION AND FUTURE WORK

We have introduced a modified formulation of the standard Warner–Terentjev energy density incorporated into a higher-order Ogden model to more accurately describe the behavior of nematic elastomers at large deformation. Furthermore, this work has initiated the study of actuation from geometries beyond flat, two-dimensional sheets by exploring a curvilinear three-dimensional geometry. We have outlined the deformation of a nematic elastomer balloon under simple expansion and twist. The material is actuated remotely by changing the temperature to dictate the degree of anisotropy, and the response is tunable. The foundation for our understanding of nematic elastomer actuation from flat geometries has already been well established with respect to the design, optimization, manufacturing, and tuning (e.g., voxelated sheets,^{28} wrinkling-resistant membranes,^{21} and moving inchworm^{31}). Future applications based on more complex geometries and loading conditions, for instance, incorporation of disclination defects and gradients of director or temperature across the thickness, can build upon this framework. Finally, the actuation can be affected by light instead of temperature.^{30} For example, the pump described here would function with a light-driven actuation from $r=8$ to $r=2$. A practical difficulty to be overcome is that many photo-active materials have a low penetration depth and the actuation is in bending, rather than stretching.

## ACKNOWLEDGMENTS

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

## DATA AVAILABILITY

The data that support the findings of this study are available within the article.

## REFERENCES

*StatPearls [Internet]*(updated May 5, 2020) (StatPearls Publishing, Treasure Island, FL, 2020).