Nematic liquid crystalline elastomers are aeolotropic materials.

Continuum models describing ideal nematic solids are widely used in theoretical studies of liquid crystal elastomers. However, experiments on nematic elastomers show a type of anisotropic response that is not predicted by the ideal models. Therefore, their description requires an additional term coupling elastic and nematic responses, to account for aeolotropic effects. In order to better understand the observed elastic response of liquid crystal elastomers, we analyse theoretically and computationally different stretch and shear deformations. We then compare the elastic moduli in the infinitesimal elastic strain limit obtained from the molecular dynamics simulations with the ones derived theoretically, and show that they are better explained by including nematic order effects within the continuum framework.

Introduction

Liquid crystalline solids are responsive multifunctional materials that combine the flexibility of polymeric networks with the nematic order of liquid crystals [1,2]. Owing to their molecular architecture, they can exhibit dramatic spontaneous deformations and phase transitions, which are reversible and repeatable under heat, light, solvents and electric or magnetic fields [3-12]. However, their physical behaviour under combined mechanical loading and external stimuli still needs to be fully elucidated.

For ideal monodomain nematic elastomers, with the mesogens uniaxially aligned throughout the material, a continuum model is given by the so-called neoclassical strain-energy density introduced in [13-16]. This is a phenomenological strain-energy function that extends the neo-Hookean model for rubber [17], where the parameters are derived from macroscopic shape changes at small strain or through statistical averaging at microscopic scale [18,19]. Since elastic stresses dominate over Frank elasticity induced by the distortion of mesogens alignment [20-22], Frank energy [23,24] is generally neglected. Extensions to polydomains where every domain has the same strain-energy density as a monodomain are provided in [25,26]. These descriptions have been generalized by employing other hyperelastic models, such as Mooney–Rivlin, Gent and Ogden, which better capture the nonlinear elastic behaviour at large strains [27-29] (see [30,31] for molecular interpretations of the Mooney–Rivlin and Gent strain energies in rubber elasticity). Further generalizations can be found in [32,33]. Numerical studies of liquid crystal elastomers (LCEs) are presented, for example, in [33-35], where the finite-element method is used, and in [36,37], where molecular Monte Carlo simulations are employed.

Usually, when the elastic properties of a material are investigated, uniaxial deformations, which are easier to reproduce experimentally, are examined first [38-41]. For a purely elastic isotropic material, the shear modulus is then inferred from a universal relation between elastic moduli from the classical theory. To study the elastic responses of nematic monodomains, uniaxial deformations were assumed in [38] (see also [42]), where it was found that, if only the elastic energy was considered, then the stretch moduli in the direction parallel to the director and in a perpendicular direction were equal. However, experiments clearly show an aeolotropic effect; namely, the stretch moduli depend on the direction in which they are measured. To capture this experimentally observed response of the material, the nematic energy [43] was then also taken into account. Experimental results for monodomains where the tensile load formed different angles with the initial nematic director were reported in [39,44]. In [45], measurements of five independent elastic constants derived from three uniaxial tests, with the director parallel, perpendicular or at an angle of 45° relative to the loading direction, respectively, were obtained for a nematic monodomain treated as a classical transversely isotropic material. However, for many complex materials, shear deformations can reveal important additional mechanical effects, which may not be observed or inferred from uniaxial tests. In particular, to assess LCE materials, shear deformations with the direction of shear either parallel or perpendicular to the nematic director need to be considered independently of uniaxial stretches [46-52]. For example, on the one hand, it was found experimentally in [38] that, if and E⊥ are the stretch moduli in a direction parallel or perpendicular to the nematic director, respectively, then at low temperature, at high temperature and at the transition point. On the other hand, if and μ⊥ denote the shear moduli in a direction parallel or perpendicular to the nematic director, respectively, then it was reported in [47,48,51] that in the isotropic phase and at temperatures below that for the nematic–isotropic phase transition. Therefore, from a symmetry point of view, monodomain LCEs are transversely isotropic materials with five independent elastic constants and the distinguished direction given by the nematic director. However, despite the constitutive symmetry about the direction given by the nematic field, the mechanical responses of LCEs differ from the known elastic behaviours in traditional transversely isotropic elastic materials where, typically, and [53-56].

The aim of this study is to develop an explicit approach for the derivation of elastic moduli that captures the aeolotropy of liquid crystalline elastomers. This approach represents an extension of the general theoretical framework by which similar elastic moduli were obtained for hyperelastic materials [57]. In the case of nematic solids, these moduli include information about both the elasticity of the polymeric network and the mechanical responses of the liquid crystal molecules. In §2, we recall the neoclassical model for ideal nematic elastomers, with the isotropic phase at high temperature as the reference configuration [29,58-61], instead of the nematic phase at cross-linking [14-16,32,33,62,63]. Phenomenologically, this choice is motivated by the multiplicative decomposition of the effective deformation into an elastic distortion, followed by a natural stress-free shape change [64-67]. This multiplicative decomposition is similar to those found in the constitutive theories of thermoelasticity, elastoplasticity and growth [68,69] (see also [70,71]), but it is also different in the sense that the elastic deformation is directly applied to the reference state. The elastic stresses can then be used to study the final deformation where the stress-free geometrical change also plays a role [65]. In §3, we calculate the two stretch moduli under small elastic uniaxial tension and finite natural deformation when the nematic director is either parallel or perpendicular to the tensile direction, respectively. In §4, we further consider three shear deformations where the elastic component is a simple shear, while the nematic director is either parallel to the shear direction, perpendicular to the shear direction or perpendicular to the shear plane, respectively. When the elastic shear strain is small and the natural deformation is finite, we obtain effective shear moduli with the relative ratio equal to the natural anisotropy parameter of the nematic material. Note that we use the uniaxial and simple shear deformation, respectively, to describe the elastic contribution to the deformation rather than the overall deformation, which also contains a natural deformation component. To derive the elastic moduli, we then take the limit of small elastic strain, while the natural deformation remains finite. This enables us to rigorously adapt the elasticity theory to nematic elastomers (for a review on elastic moduli, see [57]). To account for the physical aeolotropy of real nematic solids, in §5, we extend the continuum model by incorporating a nematic energy, and show how the stretch and shear moduli corresponding to the ideal case are modified by the additional information. The parameters entering the LCE model characterize either the change of the microstructure due to nematic effects or the behaviour of the material in large natural deformations. Nevertheless, in the limit of small elastic deformations and for fixed nematic parameters, the system behaves indeed like a transverse isotropic material with five independent elastic constants. In physics, aeolotropy refers to materials exhibiting different properties depending on the direction in which they are measured, or simply defined by Lord Kelvin as ‘That which is different in different directions’ [72, p. 122]. While traditional anisotropic elastic materials also exhibit aeolotropy, we refer to nematic solids as aeolotropic materials, and reserve the characterization of isotropic or anisotropic for the elastic part of the energy. In §6, we present a molecular dynamics simulation of a nematic elastomer, and analyse its response under similar stretch and shear deformations as for the continuum model to illustrate the aeolotropic mechanical responses. In §§3–5, physical quantities are treated symbolically, and units of measure only appear in §6, where datasets are used to illustrate the theory. The final section containsconcluding remarks.

Prerequisites

The strain-energy density describing an ideal monodomain nematic liquid crystalline (NLC) solid takes the general form [64-67,73] where F represents the deformation gradient from the isotropic state, n is a unit vector, known as the director, for the orientation of the nematic field and W(A) denotes the strain-energy density of the isotropic polymer network, depending only on the (local) elastic deformation tensor A. The tensors F and A satisfy the following relation: where is the spontaneous deformation tensor describing a change of frame of reference from the isotropic to a nematic phase (e.g. [74]). In (2.3), a > 0 represents a temperature-dependent stretch parameter, ν is the optothermal analogue to the Poisson ratio [74] and relates responses in directions parallel or perpendicular to the director n, ⊗ denotes the tensor product of two vectors and I = diag(1, 1, 1) is the identity tensor. It is assumed here that a and ν are spatially independent. The ratio r = a1/3/a− = a( represents the anisotropy parameter, which, in an ideal nematic solid, is the same in all directions. In the nematic phase, both the cases with r > 1 (prolate molecules) and r < 1 (oblate molecules) are possible, while when r = 1 the energy function reduces to that of an isotropic hyperelastic material. Nematic elastomers have ν = 1/2, while for nematic glasses ν ∈ (1/2, 2) [10,75]. Natural strains in NLC glasses are typically of up to 4%, whereas for NLC elastomers these may be up to 400%. The nematic director n is an observable (spatial) quantity. Denoting by n0 the reference orientation of the local director corresponding to the cross-linking state, n may differ from n0 both by a rotation and by a change in r. In nematic elastomers, which are weakly cross-linked, the director can rotate freely, and the material exhibits isotropic mechanical effects. In nematic glasses, which are densely cross-linked, the director n cannot rotate relative to the elastic matrix, but changes through convection due to elastic strain and satisfies [21,22,58,74,76] This constraint enables patterning of the director field at cross-linking and guarantees that the ‘written-in’ pattern remains virtually the same during natural shape changes [21,22,77].

For a hyperelastic material described by the strain-energy density W = W(A), the Cauchy stress tensor is equal to where the so-called hydrostatic pressure p denotes the Lagrange multiplier for the incompressibility constraint , β1 = 2∂W/∂I1 and β−1 = −2∂W/∂I2 are material parameters, B = AA is the left Cauchy–Green elastic deformation tensor and I1, I2 are its first two principal invariants (I3 = 1 owing to incompressibility). The corresponding first Piola–Kirchhoff stress tensor is equal to where . For an NLC solid characterized by the strain-energy density given by (2.1), the stress tensors when the director is ‘free’ to rotate relative to the elastic matrix and when the nematic director is ‘frozen’ and satisfies condition (2.4) are presented in appendix A (see also [65]).

In the next sections, we obtain two stretch moduli under small elastic uniaxial tension and three nonlinear shear moduli under elastic simple shear deformations, respectively, which combine elastic and nematic effects. In our calculations, the nematic director is ‘frozen’, but the case where the director is ‘free’ can be treated similarly, provided that the elastic deformation is small. In particular, when the director is free to rotate and a tensile force is applied perpendicular to the director, experimental results show that there is a range of strains, up to 10% (e.g. [78,79]), before the director rotates in response to the applied force. However, the local nematic order might be altered. We further note that, in finite elasticity, the strain-energy density W(A) and the stress relationship (2.5) characterize an isotropic material. Yet, the response of a nematic solid also depends on the director orientation. For instance, we show that the ideal model exhibits an anisotropic response under stretch if the natural deformation is finite. There are, in fact, two possible contributions to aeolotropy. First, in the ideal case, the elastic energy of the system is isotropic but the full energy depends on the orientation of the nematic field as well [59]. Second, there is a component of the energy that depends directly on the nematic order through the so-called Q-tensor [23]. In this case, we need to extend the discussion by including this nematic energy density as developed later in the paper.

Stretch moduli

The stretch modulus of a homogeneous isotropic elastic material is obtained under uniaxial tension with the gradient tensor in a Cartesian system of coordinates taking the form [57] where λ > 1 is the stretch ratio in the direction of the applied tensile force. Assuming that the only non-zero component of the associated first Piola–Kirchhoff stress, given by (2.6), is in the tensile direction, i.e. P22 > 0, it follows that The Young modulus at small strain is then defined as [57] where P22 is the tensile first Piola–Kirchhoff stress given by (3.2), λ − 1 is the corresponding tensile strain (for different definitions of an elastic strain, see [57]) and μ = lim (β1 − β−1) is the shear modulus at small strain.

To derive stretch moduli for the nematic material, we apply the tensile force in a direction that is either parallel or perpendicular to the reference nematic director [45]. In each case, we assume an overall deformation where the elastic component is a uniaxial tension with the deformation tensor given by (3.1). Then, we calculate the stretch moduli by taking the ratio between the first Piola–Kirchhoff tensile stress and the associated strain in the limit of small elastic tensile strain, while the natural deformation remains finite.

Nematic director parallel or perpendicular to the tensile direction

When the tensile force is acting in the second Cartesian direction and the reference director is parallel to the tensile force, the nematic director in the current configuration and the associated spontaneous deformation tensor take the following form, respectively (figure 1a): When the reference director is perpendicular to the tensile force, the nematic director in the current configuration and the associated spontaneous deformation tensor are, respectively (figure 2a), Assuming that the elastic deformation tensor A is of the form given by (3.1), for the deformation gradient F and the associated first Piola–Kirchhoff stress, the principal components in the direction of the tensile load are, respectively, and where P22 is the first Piola–Kirchhoff stress given by (3.2). These are illustrated in figures 1b and 2b, respectively, for μ = 1 and ν = 1/2. We define the following stretch moduli for the nematic material under the above two deformations, respectively: and where E = 3μ is Young’s modulus defined by (3.3). Therefore, the stretch moduli given by (3.8) and (3.9) satisfy the following relation: We infer that: if r > 1, then E(1) < E(2); if r < 1, then E(1) > E(2); if r = 1, then E(1) = E(2) = E.

Figure 1.

(a) Nematic material with the director parallel to the applied tensile force, and (b) the effect of varying the parameter a on the first Piola–Kirchhoff tensile stress of an ideal material when μ = 1 and ν = 1/2. (Online version in colour.)

Figure 2

(a) Nematic material with the director perpendicular to the applied tensile force, and (b) the effect of varying the parameter a on the first Piola–Kirchhoff tensile stress of an ideal material when μ = 1 and ν = 1/2. (Online version in colour.)

Shear moduli

To obtain the shear modulus of a homogeneous isotropic elastic material, a standard deformation is the simple shear, with the following gradient tensor in a Cartesian system ofcoordinates [57]: where k > 0 is the shear parameter. The non-zero components of the associated Cauchy stress tensor, given by (2.5), are and the non-zero components of the corresponding first Piola–Kirchhoff stress tensor, given by (2.6), are The shear modulus at small shear is then obtained as follows [57]: To obtain suitable shear moduli for the nematic solid, we assume an overall deformation where the elastic component is a simple shear with the deformation tensor given by (4.1). In general, one cannot easily separate the individual contributions to the deformation gradient given by (2.2). However, for the three shear configurations adopted here, it is possible to write these components explicitly. Then, in each case, we can calculate the shear modulus by taking the ratio between the first Piola–Kirchhoff shear stress and the corresponding shear strain in the limit of small elastic shearing strains, while the natural deformation remains finite.

Nematic director parallel to the shear direction

When the reference nematic director is parallel to the direction of applied shear force, in the current configuration, we have (figure 3a) By (2.2) and (A4), the corresponding shear strain and first Piola–Kirchhoff shear stress for the nematic solid are, respectively, where P12 is the shear component of the elastic first Piola–Kirchhoff stress given by (4.3). These are represented in figure 3b, for μ = 1 and ν = 1/2. We now define the shear modulus for the nematic material at small shear as follows:

Figure 3

(a) Nematic material with the director parallel to the applied shear force, and (b) the effect of varying the parameter a on the first Piola–Kirchhoff shear stress when μ = 1 and ν = 1/2. (Online version in colour.)

Nematic director perpendicular to the shear direction

When the reference nematic director is perpendicular to the direction of shear, we obtain in the current configuration (figure 4a) For sufficiently small values of k, the corresponding shear strain and first Piola–Kirchhoff shear stress of the nematic material take the form with the components P12 and P22 of the elastic first Piola–Kirchhoff stress given by (4.3). These are illustrated in figure 4b, for μ = 1 and ν = 1/2. In this case, the associated shear modulus at small shear is equal to

Figure 4

(a) Nematic material with the director perpendicular to the applied shear force, and (b) the effect of varying the parameter a on the first Piola–Kirchhoff shear stress when μ = 1 and ν = 1/2. (Online version in colour.)

Nematic director perpendicular to the shear plane

When the reference nematic director is perpendicular to the shear plane, we have (figure 5a) By (2.2) and (A4), the corresponding shear strain and first Piola–Kirchhoff shear stress for the nematic solid are, respectively, where P12 is the shear component of the elastic first Piola–Kirchhoff stress given by (4.3). These are represented in figure 5b, for μ = 1 and ν = 1/2. We now define the shear modulus for the nematic material at small shear as follows: A comparison of the shear moduli given by (4.7), (4.10) and (4.13) implies We conclude that: if r > 1, then μ(1) < μ(2) < μ(3); if r < 1, then μ(1) > μ(2) > μ(3); if r = 1, then μ(1) = μ(2) = μ(3) = μ.

Figure 5

(a) Nematic material with the director perpendicular to the shear plane, and (b) the effect of varying the parameter a on the first Piola–Kirchhoff shear stress when μ = 1 and ν = 1/2. (Online version in colour.)

Contribution of the nematic free energy

The results of §§3 and 4 imply that, for an ideal nematic material, the effective shear and stretch moduli respect the same inequalities (e.g. if r > 1, then we have both E(1) < E(2) and μ(1) < μ(2) < μ(3), and so on). However, numerous experimental results have demonstrated that there are significant differences between the behaviour of real nematic solids and that of ideal nematic materials analysed in the previous sections [38,42,45]. In particular, it was found that the ideal behaviour of the stretch and shear moduli ratio does not match experimental results. We therefore extend the strain-energy function by taking into account the nematic free-energy density, in addition to the isotropic strain-energy density W(nc) = W(nc)(F, n), given by (2.1), as where W is equal to [38,43] (see also [19], ch. 2) with Q the order parameter tensor. This macroscopic tensor parameter is used to describe the orientational order in nematic liquid crystals [23] (see also [80]). For incompressible nematic elastomers subjected to uniaxial stretches, the contribution given by (5.2) to the total strain-energy density described by (5.1) was originally analysed in [38]. Following a similar approach, we restrict our attention to the one-term strain-energy density of the form (2.1) (see also [58]), where ‘tr’ denotes the trace operator and μ > 0 is the elastic shear modulus at small strain. This strain-energy density can be expressed equivalently as follows [14]: where is the natural deformation tensor corresponding to the cross-linking state (with eigenvalues and , denoting the direction parallel to the director and ⊥ indicating an orthogonal direction), is the natural deformation tensor in the current configuration (with eigenvalues , , ) and satisfies the relation [29]. When the order parameter tensor takes the form Q = diag( − (Q − b)/2, − (Q + b)/2, Q), where Q and b are scalar values (see also [19], §2.2), we consider the nematic energy described by (5.2) and approximate it by [38] We therefore approximate the total strain-energy density defined by (5.1) as where and are given by (5.4) and (5.5), respectively. For the components of G2, the first-order approximation of the Taylor expansions about the backbone order parameters in the initial state (Q, b) = (Q0, 0) are [38] where and l⊥, = ∂l2/∂Q = ∂l3/∂Q denote the first derivatives of l1 and l2 (and also l3) with respect to Q at b = 0, respectively, and l2, = ∂l2/∂b and l3, = ∂l3/∂b are the first derivatives of l2 and l3, respectively, with respect to b at Q = Q0.

For a nematic solid with the strain-energy density given by (5.6), we derive the stretch and shear moduli in a direction parallel or perpendicular to the nematic director, as follows.

Extension parallel to the director

When the nematic director is uniformly aligned in the first direction, we take , where λ1 = λ > 1 and . Assuming infinitesimal extension, we have λ = 1 + ε and λ−1 = 1 − ε, where ε denotes the infinitesimal strain. At b = 0, the elastic strain-energy density is approximated by This contributes to the small variation in total strain-energy density defined by (5.6) where . The above quadratic function in δQ has a minimum of Denoting E = 3μ, the stretch modulus in a direction parallel to the nematic director is then

Extension perpendicular to the director

When the director is aligned in the first direction and , where and λ2 = λ > 1, assuming λ = 1 + ε and λ−1 = 1 − ε, with ε the infinitesimal strain, the elastic strain-energy density given by (5.4) is approximated by This contributes to the total strain-energy variation where is defined as before and . Minimizing the above function with respect to δQ and b gives Denoting E = 3μ, the stretch modulus in a direction perpendicular to the nematic director is

Shear parallel to the director

We recall that a simple shear deformation of a hyperelastic material is equivalent to a biaxial stretch (‘pure shear’ [81]) in the principal directions [82]. We assume that the director is aligned uniformly in the first (or second direction), and , where λ1 = λ > 1, λ2 = 1/λ and λ3 = 1, and R is the rotation by an angle of π/4 in the plane formed by the first two directions. Taking λ = 1 + ε and λ−1 = 1 − ε, with ε the infinitesimal strain, the elastic strain-energy density is approximated as follows: This contributes to the total strain-energy variation with and defined as before. We minimize this function with respect to δQ and b, The linear shear modulus in a direction parallel to the nematic director is then

Shear perpendicular to the director

Next, assuming that the director is aligned in the first direction, we set , where λ1 = 1, λ2 = λ > 1 and λ3 = 1/λ. Taking λ = 1 + ε and λ−1 = 1 − ε, with ε the infinitesimal strain, the elastic strain-energy density is approximated by This contributes to the small variation in total strain-energy density where is defined as before. The above quadratic function in b has a minimum value of The linear shear modulus in a direction perpendicular to the nematic director is equal to

Deviation from isotropy

The stretch and shear moduli , E⊥, and μ⊥, defined by (5.11), (5.15), (5.19) and (5.23), respectively, contain information about the nematic order, in addition to the small strain elasticity of the polymer matrix, described by Young’s modulus E and shear modulus μ. Proceeding as in §§3 and 4, with instead of W(nc), by replacing E with in (3.8) and with E⊥ in (3.9), μ with in (4.7) and (4.10), and with μ⊥ in (4.13), we obtain We deduce that: if , then E(2)/E(1) < r2; if , then E(2)/E(1) > r2; if , then μ(3)/μ(2) < r; if , then μ(3)/μ(2) > r. In particular, when r > 1, if and , then E(1) > E(2) and μ(1) < μ(2) < μ(3), which is qualitatively different from the behaviour of ideal nematic solids or any standard anisotropic hyperelastic material. In the next section, we show how to access these moduli through molecular dynamics simulations.

Molecular dynamics simulation

We performed molecular dynamics (MD) simulations to synthesize LCEs, form nematic LCEs and characterize their response under stretch and shear deformations. Given its accuracy and efficiency for modelling mesogen–polymer systems, we used a hybrid force field, including Gay–Berne coarse-graining potentials for mesogen–mesogen interaction and Lennard-Jones (LJ) potentials for united atoms of hydrocarbon groups, CH, in polymer chains. Our computer simulations [83] can serve as a virtual experiment to observe the macroscopic mechanical behaviour, which can then be compared directly with the continuum theory, and to calculate the physical quantities arising from atomistic movement and gain a mechanistic understanding.

In the MD simulations of LCEs, the potential energy of the whole system includes contributions from bond stretch, angle bending, dihedral rotation, non-bonded LJ interaction between united atoms (a-a), anisotropic non-bonded Gay–Berne interaction between mesogens (m-m) and the extended Gay–Berne interaction between united atoms and mesogens (a-m): where In the above equations, Ebond represents a harmonic bond style, with l the ith bond length, l(1) the equilibrium bond length of CH–CH and l(2) the equilibrium bond length of CH–mesogen; Eangle denotes a harmonic angle style, with θ the ith angle, θ(1) the equilibrium angle of the non-branched X–CH2–X and θ(2) the equilibrium angle of branched X–CH–X; Edihedral is a multiple-term harmonic dihedral style, with ϕ the ith torsion angle, the non-branched X–CH2–CH2–X and the branched X–CH2–CH–X; E( represents the LJ potential between non-bonded united atoms CH, with r the distance between two united atoms and a, c the factorized energy parameters for CH; E is the Gay–Berne potential between non-bonded mesogens, with U the shifted distance-dependent interaction, η, χ the orientation-dependent energy, A the transformation matrix for mesogen i, r the centre-to-centre vector between the ith and jth mesogens and all the rest of the parameters specified as constants in table 1; and E denotes the extended Gay–Berne potential between a non-bonded mesogen and CH following the standard mixing rule [86]. The cut-off distance is 9.8 Å for the LJ potential and 16.8 Å for the Gay–Berne potential. For the detailed explanation and explicit form of the Gay–Berne potential, readers should refer to [87,88]. First, 64 molecules of side-chain liquid crystal polymers were created, where every molecule has a backbone of 100 hydrocarbon monomers and 50 side chains attaching to the backbone in a syndiotactic way. Among the 50 side chains for each molecule, 20% were randomly selected to be attached with cross-linking sites, and the rest were attached with mesogens, as shown in figure 6a. Thus, every LC molecule has a different configuration. The 64 molecules were mixed by heating at 800 K for 50 ns. Then the system was quenched down to 500 K during 10 ns. Cross-linking of the first step was established while equilibrating the system at 500 K. A weakly cross-linked isotropic LCE was constructed, as shown in figure 6b. The system was found to have its isotropic–nematic phase-transition temperature below 490 K, consistent with other MD studies of similar LCE systems [89]. To form nematic LCEs, the whole system was quenched down to 450 K with an external field , where P2(x) = (3x2 − 1)/2 is the second Legendre polynomial and θ is the angle between the long axis of the ith mesogen and the external field. The external field has been experimentally and computationally proved to accelerate the formation of the nematic phase [90,91]. The quenching stage from 500 K to 450 K lasted 10 ns and the equilibrium stage at 450 K lasted 20 ns, for both of which the external field was applied along the z-direction. The external field was removed while the LCE system was equilibrated at 450 K for another 17 ns and the nematic phase was found to be stable. The nematic order S2 = 〈P2(cosθ)〉 is shown in figure 6b during the quenching and equilibrium state. Subsequently, the second stage of cross-linking was performed within the nematic LCEs. The nematic system was then quenched down to 300 K and the nematic order was found to be well maintained. Throughout the whole process of constructing nematic LCEs, NPT calculations were performed using a time step of 1 fs in the Large-scale Atomic/Molecular Massively Parallel Simulator (LAMMPS) with periodic boundary conditions along three dimensions. At the isotropic–nematic transition, a spontaneous deformation was observed in the MD simulations shown in figure 6b, owing to the alignment of ellipsoidal mesogens along the Z-direction. The nematic LCE at 300 K has three dimensions: , and . During the isotropic–nematic phase transition, the volume of LCEs should demonstrate a negligible change [92]. However, the MD simulations of LCEs using Gay–Berne potentials are known to show thermal expansion effects [93]. To eliminate these non-physical effects, the isotropic LCE at 500 K was quickly quenched down to 300 K and the three dimensions for the isotropic phase were measured as , and . The respective stretch ratios were then calculated as λ = 0.8280, λ = 0.7426 and λ = 1.5844. Ideally, the biaxial Gay–Berne potential for mesogen–mesogen interaction should yield the same stretch ratio along the X- and Y-directions. Here, however, the small size of the LCE system in the MD simulations gives rise to slightly different λ and λ. To eliminate the size effect on the spontaneous deformation, these stretch ratios were averaged to λ = 0.7853. Then the anisotropy ratio can be estimated as r ≈ λ/λ = 2.0176, Poisson’s ratio is ν ≈ 0.5 and a = r2 ≈ 4.0706.

Table 1.

Parameters for interatomic potentials used in the MD simulation [84,85].

parameters	value (units)
kb, bond energy constant	520.0156 (kcal/mol/Å2)
l(1), equilibrium CHx–CHy bond length	1.540 (Å)
l(2), equilibrium CHx–mesogen bond length	7.075 (Å)
ka, angle energy constant	124.2009 (kcal/mol)
θ(1), equilibrium non-branched angle	114.0014 (°)
θ(2), equilibrium branched angle	112.0018 (°)
C0(1), non-branched torsion energy constant	2.0066 (kcal/mol)
C0(2), non-branched torsion energy constant	4.0111 (kcal/mol)
C0(3), non-branched torsion energy constant	0.2709 (kcal/mol)
C0(4), non-branched torsion energy constant	−6.2885 (kcal/mol)
C1(1), branched torsion energy constant	0.7413 (kcal/mol)
C1(2), branched torsion energy constant	1.8264 (kcal/mol)
C1(3), branched torsion energy constant	0.5329 (kcal/mol)
C1(4), branched torsion energy constant	−3.4521 (kcal/mol)
ai, energy parameter for CH3	2534.0341 (kcal/mol⋅Å^6)
ci, energy parameter for CH3	47.2921 (kcal/mol⋅Å^3)
ai, energy parameter for CH2	2251.9322 (kcal/mol⋅Å^6)
ci, energy parameter for CH2	37.0955 (kcal/mol⋅Å^3)
ai, energy parameter for CH1	1467.0522 (kcal/mol⋅Å^6)
ci, energy parameter for CH1	21.2860 (kcal/mol⋅Å^3)
ϵ, well depth for Ur function in (6.1a–f)	0.8079 (kcal/mol)
σ, minimum effective particle radii in (6.1a–f)	5 (Å)
ϵi, relative well depth for end-to-end	5
ϵi, relative well depth for side-to-side	1
v, exponent for η function in (6.1a–f)	1
ξ, exponent for χ function in (6.1a–f)	2
κ, length/breadth ratio of mesogen	3
mass, united atoms CHx	12.0 + x (g/mol)
mass, mesogen	226.0 (g/mol)

Figure 6.

(a) Schematic of single chains of a liquid crystal polymer with a hydrocarbon backbone and 50 side chains among which 20% are attached with cross-linking sites (yellow atoms) and 80% are attached with mesogens (white ellipsoids). The cross-linking sites and mesogens are randomly selected from the 50 side chains for each liquid crystal polymer molecule and chain 1 and chain 2 are displayed here as examples. (b) The evolution of nematic order S2 during the isotropic–nematic phase transition when the weakly cross-linked LCEs were first quenched from 500 K to 450 K with external field for 10 ns, then equilibrated at 450 K with external field for 20 ns, then equilibrated at 450 K without external field for 17 ns. (Online version in colour.)

After obtaining the nematic LCE at 300 K, we subjected the system to either stretch or shear deformations as described in §§3 and 4, respectively. For each deformation, we can calculate the evolution of the first Piola–Kirchhoff stresses with respect to the strains. We first calculate the Cauchy stress tensor in LAMMPS based on , where m represents the mass of the kth atom, v denotes the velocity of the kth atom in the ith dimension, N is the total number of atoms, V is the total volume of the system, N′ is the number of atoms pairs and r and f are the positions and forces of atom k along the ith and jth dimensions, respectively [94]. Next, given the deformation imposed, we construct the elastic deformation tensor A and, subsequently, the corresponding first Piola–Kirchhoff stress tensor P, following (2.6). Then, imposing the spontaneous deformation tensor G, or simply following (3.6), (3.7), (4.6), (4.9) or (4.12), we obtain the final . The results are shown in figures 7 and 8, respectively, where all the cases demonstrate some nonlinearity in the stress-deformation relations. From the two stretch deformations, we derive the effective stretch moduli E(1) = 2253.4 MPa and E(2) = 775.475 MPa. From the three shear deformations, we find the effective shear moduli μ(1) = 132.78 MPa, μ(2) = 258.26 MPa and μ(3) = 889.97 MPa. Thus, These relations are in agreement with (5.24) when and (see also fig. 7 of [38], fig. 7 of [45], fig. 6 of [48] and fig. 3 of [51]). Our theoretical results are independent of the manufacturing history. Experimentally, the mechanical properties of main-chain and side-chain LCEs are compared in [51]. In future work, it would be interesting to carry out separate simulations for these different cases to compare their behaviour computationally as well.

Figure 7.

The first Piola–Kirchhoff tensile stress in the nematic LCE system when (a) the director is parallel to the tensile force and (b) the director is perpendicular to the tensile force. (Online version in colour.)

Figure 8.

The first Piola–Kirchhoff shear stress in the nematic LCE system when (a) the director is parallel to the shear force, (b)the director is perpendicular to the shear force and (c) the director is perpendicular to the shear plane. (Online version in colour.)

Conclusion

We studied theoretically and computationally the mechanical behaviour of nematic LCEs under different stretch and shear deformations. Theoretically, we first examined ideal nematic elastomers characterized by a homogeneous isotropic elastic strain-energy density, then also phenomenological models incorporating an additional nematic energy. We showed that these cases are qualitatively different, and that the generalized model does not necessarily order stretch moduli in the same way as the shear moduli. We also performed molecular dynamics simulations to analyse numerically the responses of simulated systems under a similar set of deformations, and found different mechanical responses in different directions. Therefore, the trifecta of experiments, computations and theory leads us to conclude that the contribution of the nematic free energy cannot be ignored, even in small deformations, and that LCEs are best understood as aeolotropic materials. When Frank effects also play an important role, they need to be taken into account as well. However, the deviation from isotropy is well captured by including the nematic energy, and this constitutes an important step in the constitutive modelling of liquid crystalline solids.