A predictive theoretical model for stretch-induced instabilities in liquid crystal elastomers

ABSTRACT Following the experimental lead, we construct a general mathematical model which, depending on whether the uniaxial scalar order parameter is constant or not, can predict either the classical shear striping instability or the molecular auxetic response and mechanical Fréedericksz transition observed in different liquid crystal elastomers. Our theoretical model can shed some light on the role of nematic order in these stretch-induced mechanical instabilities not observed in conventional rubber-like solids. Graphical abstract


Introduction
In Refs [1][2][3][4][5][6], experimental observations are reported for a nematic liquid crystal elastomer (LCE) where an auxetic response occurs, such that, when a material sample is stretched longitudinally, its volume remains unchanged, but its thickness, measured in the direction perpendicular to the sample's plane, first decreases, then increases if the deformation is sufficiently large.This behaviour is accompanied by a mechanical Fréedericksz transition whereby the nematic director, which is initially aligned within the sample's plane and perpendicular to the longitudinal direction, is mechanically deformed to become parallel to the applied force.Experimentally, this can appear as a sudden rotation of the director, and is different from the gradual rotation occurring in many other LCEs where alternating shear stripes form at very low stress [7][8][9][10][11][12][13]. Discontinuous rotation of the nematic director was also reported previously by Mitchell et al. [14] (see [15,Section 5.5]).
In Mihai et al. [25], a descriptive three-term Ogdentype model [26,27] is adopted which permits the calibration to available experimental data for the auxetic LCE.However, despite the inherent mathematical simplicity of the Ogden strain-energy function, fewer terms are needed in order to assess the effect of changing parameter values.
In this paper, guided by experimental observations, we propose a two-term Ogden-type constitutive model, introduced in Section 2, which, depending on whether the uniaxial scalar order parameter is constant or not, predicts either the auxetic response and the mechanical Fréedericksz transition (Section 2.1) or the classical shear striping (Section 2.2) presented by various LCEs (see Figure 1).These results are illustrated by a set of examples presented in Section 3.

Model function
We consider the following two-term strain-energy function describing a nematic LCE (see also [25]), where μ 1 > 0, μ 2 > 0 and n > 0 are constant parameters, and μ ¼ μ 1 þ μ 2 represents the shear modulus at infinitesimal strain.In the above equation, fλ 2 1 ; λ 2 2 ; λ 2 3 g denote the eigenvalues of the Cauchy-Green tensor FF T , where F is the deformation gradient from the reference crosslinking state, and fα 2 1 ; α 2 2 ; α 2 3 g are the eigenvalues of AA T , where A ¼ G À 1 FG 0 is the local elastic deformation tensor, with G 0 and G the 'natural' deformation tensors due to the liquid crystal director in the reference and current configuration, respectively.These tensors satisfy the following relations [29] (see also [15,Chapter 3]): with c 0 and c representing the effective step length of the polymeric chain, Q 0 and Q the symmetric traceless order parameter tensors [15, pp.48-49], and I ¼ diagð1; 1; 1Þ the tensor identity.The macroscopic tensor parameters are known to describe orientational order in nematic liquid crystals [30].
In the following sections, we demonstrate that, for a suitable choice of the constitutive parameters, the model defined by ( 1) is capable of predicting either a large-strain auxetic response accompanied by mechanical Fréedericksz transition or the shearstriping instability.

Variable order parameter and mechanical Fréedericksz transition
First, we examine the emergence of large-strain auxetic behaviour.To obtain this, for the LCE sample, in a Cartesian system of coordinates ðX 1 ; X 2 ; X 3 Þ, we designate the plane formed by the first two directions as the sample's plane and the third direction as its thickness direction.The material sample is elastically deformed by application of a tensile force in the first (longitudinal) direction, while the nematic director is initially along the second (transverse) direction.
Setting the nematic director in the reference and current configuration, respectively, as follows: where θ 2 ½0; π=2� is the angle between n and n 0 , the deformation gradient takes the form with λ 1 λ 2 λ 3 ¼ 1 and diagð�; �; �Þ denoting a diagonal second-order tensor.
Guided by experimental observations (see also [24]), we assume that the rotation of the nematic director within the plane formed by its initial orientation and the applied force from θ ¼ 0 to θ ¼ π=2 happens suddenly [25], and denoted by λ crt the longitudinal stretch ratio λ where the rotation occurs.In practice, although the angle θ will also take intermediate values between 0 and π=2, the deformation interval for which the director rotation is observed is very short, and its separate analysis can be omitted.However, in this case, a jump may be observed in the orthogonal stretch ratios λ 2 and λ 3 at the critical longitudinal stretch λ 1 ¼ λ crt .
In the reference configuration, the LCE is uniaxial, and the order parameter tensor is equal to [15, p. 14] where Q 0 is the scalar order parameter.
In the deformed configuration, biaxiality can also emerge [15, p. 15], and therefore • If θ ¼ 0, then the order parameter tensor takes the form where Q and b are the uniaxial and biaxial scalar order parameters, respectively.Numerically, we will consider both the cases when biaxiality is neglected, i.e. b ¼ 0, and when the magnitude of biaxial parameter b increases (decreases) as the uniaxial order parameter decreases (increases).
For the corresponding elastic Cauchy-Green tensor , where ''T 00 denotes the transpose (see also [25]), we have Under the uniaxial deformation considered here, the second and third directions must be stress free.Then the corresponding principal Cauchy stresses, defined by with p denoting the usual Lagrange multiplier for the incompressibility constraint, take the form The associated first Piola-Kirchhoff stresses are Taking λ 1 ¼ λ, we solve for λ 2 the equation and obtain To capture the auxetic response, while keeping our continuum mechanics formulation tractable, we let the uniaxial scalar parameter Q > 0 decrease linearly with respect to the longitudinal stretch ratio λ from a given value Q 0 2 ð0; 1Þ, at λ ¼ 1, to a minimum value Q crt 2 ð0; Q 0 Þ that occurs at the critical stretch ratio λ ¼ λ crt , then increase again, also linearly, to a value Similarly, we let the biaxial scalar parameter b < 0 to decrease linearly with respect to the longitudinal stretch ratio λ from a given value b 0 2 ðÀ 1; 0Þ, at λ ¼ 1, to a minimum value b crt 2 ðÀ 1; b 1 Þ, occurring also at the critical stretch ratio λ ¼ λ crt , then increase linearly to a value b Note that the magnitude (absolute value) of this parameter increases then decreases while its numerical value decreases then increases.
To find the critical stretch ratio λ ¼ λ crt , we solve the equation Once λ crt is obtained, we can also determine the minimum stretch ratio λ ¼ λ aux at which the auxetic behaviour begins to be observed, i.e.where λ 3 starts to increase.
The purpose of this paper is to investigate numerically how a theoretical model of the form given by Equation (1) can capture either a large strain auxetic effect similar to that shown in Mihai et al. [25] or the classical shear striping.After several numerical tests, the proposed model with n ¼ 2 was selected as being sufficiently simple to analyse in some detail and able to present each of those effects under different conditions, as demonstrated in the following sections.Henceforth, we set n ¼ 2 in the model strain-energy function.

Constant order parameter and shear striping
Next, we analyse the deformation of the LCE model defined by Equation (1) under uniaxial tensile force, such that the nematic director is uniformly distributed in the sample's plane and rotates slowly from the initial position perpendicular to the tensile force to align with this force via an intermediate state where alternating shear stripes occur.For the corresponding deformation, the stretch ratios in the directions perpendicular to the applied force are equal, agreeing with experimental observations recorded, for example, in Figure 6 of [31] and Figure 6 of [32].
In this case, we assume that the scalar order parameter Q remains constant and there is no biaxiality, i.e. b ¼ 0. This assumption is based on experimental observations showing that variations in the order parameter are relatively small before, during and after stripe domains form [8,12,13,32].The natural gradient tensor then takes the form where represents the natural anisotropy parameter for the nematic solid, n is the nematic director, � denotes the tensor product of two vectors, and I is the identity tensor.
In particular, when the director for the reference and current configuration are given by (3), the associated natural deformation tensors described by (2) are, respectively, and To demonstrate shear-striping instability, we consider the following perturbed deformation gradient where λ is the stretch ratio in the longitudinal direction, which is the direction of the applied tensile force, and 0 < ε � 1 is a small shear parameter.
The elastic deformation tensor A ¼ G À 1 FG 0 is then equal to Taking n ¼ 2 in the strain-energy function described by Equation ( 1), we obtain where the eigenvalues fλ 2 1 ; λ 2 2 ; λ 2 3 g of the Cauchy-Green tensor FF T and fα 2 1 ; α 2 2 ; α 2 3 g of the tensor AA T satisfy the following relations, respectively, and with "tr"denoting the trace and "Cof"the cofactor.Hence, and We further define the following function: with W ðlceÞ ðλ 1 ; λ 2 ; λ 3 ; Q; b; θÞ ¼ W ðlceÞ described by Equation (33).Differentiating the above function with respect to ε and θ, respectively, gives and As the equilibrium solution minimises the energy, it must satisfy the simultaneous equations At ε ¼ 0 and θ ¼ 0, both the partial derivatives defined by Equations ( 39) and (40) are equal to zero, i.e.Therefore, this trivial solution is always an equilibrium state, and for sufficiently small values of ε and θ, we have the second-order approximation where First, we find the equilibrium value θ 0 for θ as a function of ε by solving the second equation in (41).By the approximation (Equation (43)) the respective equation takes the form and implies Next, substituting θ ¼ θ 0 ðεÞ in (43) gives the following function of ε, Depending on whether the expression on the right-hand side in Equation ( 49) is positive, zero, or negative, the respective equilibrium state is stable, neutrally stable or unstable [23].
We deduce that the equilibrium state with ε ¼ 0 and θ ¼ 0 is unstable if where λ � solves the equation Similarly, at ε ¼ 0 and θ ¼ π=2, both the partial derivatives defined by Equations ( 39) and (40) are equal to zero, and Thus the equilibrium state with ε ¼ 0 and θ ¼ π=2 is unstable if where λ � solves the equation There is also an inhomogeneous equilibrium solution ðε 0 ; θ 0 Þ when λ satisfies either (50) or (55).For the equilibrium angle θ 0 ðεÞ, given by Equation (48), the second equation in (41) is satisfied.To find the equilibrium value ε 0 of ε as a function of λ, we solve the simultaneous equations Equation (41) numerically.

Results
To illustrate the mechanical behaviour of the model described by Equation (1), first we assume that the nematic director rotates suddenly, as described in Section 2.1, and consider the following four cases: (I) Variable uniaxial order parameter Q and variable biaxiality parameter b, given by Equations ( 23) and (24), respectively; (II) Variable uniaxial order parameter Q given by Equation ( 23), without biaxiality (b ¼ 0); (III) Constant uniaxial order parameter Q and variable biaxiality parameter b given by Equation ( 24); (IV) Constant uniaxial order parameter Q, without biaxiality (b ¼ 0).
We then turn our attention to case where the director rotates slowly and the uniaxial order parameter Q is constant, as analysed in Section 2.2.
In our numerical examples, we chose the maxima and minima of the order parameters Q and b similar to those reported in previous papers where experimental results were described.From those values, the coefficients of the respective piecewise linear functions are inferred, as follows: In each case, we set the same initial model parameters of μ 1 ¼ 5 μ 2 ¼ 1, n ¼ 2, and Q 0 ¼ 0:5880, the latter being a typical value for the scalar uniaxial order parameter in a nematic liquid crystalline material.
Figure 5 shows the corresponding transverse stretch ratios λ 2 and λ 3 as functions of the longitudinal stretch ratio λ 1 ¼ λ.This figure provides insight into the nature of the auxetic response seen in some LCEs.In case (I), we see how non-constant Q and b lead to a physically realistic deformation.Critically, this case predicts large anisotropy in the transverse stretch ratios λ 2 and λ 3 , and also a physically realisable auxetic response.In case (II), where Q varies and b ¼ 0, a similar behaviour is predicted, but with a less pronounced auxetic response, as the interval between λ aux and λ crt is smaller than in case (I).In cases (III) and (IV), although an auxetic reponse is predicted at λ crt , the model shows large discontinuity in λ 2 and λ 3 at The sine function of the angle θ between the nematic director in the current and reference configuration for cases (I)-(IV), respectively.For cases (I) and (II), the vertical lines correspond to the predicted longitudinal stretch ratio λ crt where the director rotates suddenly and the minimum stretch ratio λ aux < λ crt where auxeticity is observed.For cases (III) and (IV), the vertical line corresponds to λ crt .λ 1 ¼ λ crt , which can be indicative of the fact that the director may rotate slowly rather than suddenly.We can also see that, before λ crt is reached, in case (III), the non-zero biaxial order parameter allows the evolution of anisotropy in the transverse stretch ratios, while in case (IV), where the uniaxial order parameters is constant, an isotropic deformation in the transverse directions is predicted.Figure 6 displays the associated first Piola-Kirchhoff stresses.
In Figures 2-6, the predicted longitudinal stretch ratio λ crt where the director rotates suddenly is indicated in each case.For cases (I) and (II), the minimum stretch ratio λ aux < λ crt where auxeticity is obtained is also presented.
For the case when the director rotates slowly, in Figure 8, we show the constant uniaxial and biaxial order parameters Q and b, respectively, and the stretch ratios λ 1 ¼ λ and λ 2 ¼ λ 3 ¼ 1= ffi ffi ffi λ p . The constant anisotropy parameter given by Equation ( 28) when Q ¼ 0:5880 is a ¼ 5:2816.

Conclusion
LCE stretching ultimately causes a reorientation of the nematic director which tends to become parallel to the applied force.Sometimes a change in order also occurs.This correlation between macroscopic deformation and molecular architecture causes a variety of mechanical instabilities not observed in conventional rubber-like solids.
For example, experimental evidence demonstrates that, for LCEs with planar mesogen alignment, when a uniaxial tensile force is imposed perpendicular to the nematic director, the director rotates either slowly giving rise to a shear stripe pattern or almost suddenly causing auxeticity at a characteristic large strain.Experimentally, for auxetic LCEs, biaxiality has also been found to emerge.
The present theoretical study proposes a general mathematical model capable of predicting either the classical shear-striping instability or the auxetic response observed in various LCEs.From the modelling point of view, these two phenomena are captured by letting the uniaxial scalar order parameter decrease continuously to a sufficiently low value then increase again in the case of auxetic LCEs, and by maintaining it constant for LCEs exhibiting shear stripes.The theoretical model further suggests that biaxiality also plays a role in enhancing the auxetic reponse.These results may shed some light on the role of nematic order in the observed mechanical behaviours.
Further experimental investigation where nematic parameters must be varied systematically is required to fully understand the auxetic phenomenon and will be a subject for our future study.

Disclosure statement
No potential conflict of interest was reported by the authors.

Funding
We are grateful for the support by the Engineering and Physical Sciences Research Council of Great Britain under research grants EP/R020205/1 to AG, EP/S028870/1 to LAM and EP/M009521/1 to HFG.DM thanks the Leverhulme Trust for an Early Career Fellowship.

Figure 2 .
Figure 2. (Colour online) Scalar uniaxial and biaxial order parameters for cases (I)-(IV), respectively.For cases (I) and (II), the vertical lines correspond to the predicted longitudinal stretch ratio λ crt where the director rotates suddenly and the minimum stretch ratio λ aux < λ crt where auxeticity is obtained.For cases (III) and (IV), the vertical line corresponds to λ crt .

Figure 3 .
Figure 3. (Colour online) The sine function of the angle θ between the nematic director in the current and reference configuration for cases (I)-(IV), respectively.For cases (I) and (II), the vertical lines correspond to the predicted longitudinal stretch ratio λ crt where the director rotates suddenly and the minimum stretch ratio λ aux < λ crt where auxeticity is observed.For cases (III) and (IV), the vertical line corresponds to λ crt .

Figure 4 .
Figure 4. (Colour online) Energy functions for cases (I)-(IV), respectively.For cases (I) and (II), the vertical lines correspond to the predicted longitudinal stretch ratio λ crt where the director rotates suddenly and the minimum stretch ratio λ aux < λ crt where auxeticity is obtained.For cases (III) and (IV), the vertical line corresponds to λ crt .

Figure 5 .
Figure 5. (Colour online) The transverse stretch ratios for cases (I)-(IV), respectively.For cases (I) and (II), the vertical lines correspond to the predicted longitudinal stretch ratio λ crt where the director rotates suddenly and the minimum stretch ratio λ aux < λ crt where auxeticity is obtained.For cases (III) and (IV), the vertical line corresponds to λ crt .

Figure 6 .
Figure 6.(Colour online) The first Piola-Kirchhoff stresses for cases (I)-(IV), respectively.For cases (I) and (II), the vertical lines correspond to the predicted longitudinal stretch ratio λ crt where the director rotates suddenly and the minimum stretch ratio λ aux < λ crt where auxeticity is obtained.For cases (III) and (IV), the vertical line corresponds to λ crt .