A theoretical liquid crystal elastomer model that mimics the elasticity of cat skin

A mathematical model for nematic liquid crystal elastomers is proposed that mimics the elastic response of cat skin where reorientation of dermal fibres produces an increase in the thickness direction under tensile stretch. To capture this unusual effect, the uniaxial order parameter in the nematic elastomer model is allowed to decrease then increase again, and the critical stretch at which this change of monotonicity occurs and where the director also rotates suddenly is predicted. In addition, the model parameters are described by probability density functions and their uncertainty is propagated numerically to the predicted mechanical results.


Introduction
Extensive experimental programs have demonstrated the capability of a class of nematic liquid crystal elastomesr (LCEs) to present certain auxetic responses, while their volume is conserved under large elastic strain [8][9][10][11][12][13]. In particular, when a material sample is stretched longitudinally, its thickness, measured in the direction perpendicular to the sample's plane, first decreases, then increases if the deformation is sufficiently large, then decreases again. Such mechanical behaviour is generated at a molecular level, without porosity emerging at macroscopic scale, and it 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, rotates suddenly to become parallel with the applied force.
Auxetic materials that are simple to fabricate and avoid porosity-related weakening have long been sought for applications in a range of industries including sports equipment, aerospace, biomedical materials and architecture. Because of the intrinsic similarities between LCEs and conventional rubber [17], suitable descriptions can be achieved by adapting existing hyperelastic models for rubber-like solids [5]. Hyperelastic materials are the class of material models described by a strain-energy density function with respect to a fixed reference configuration.
In [6], a quantitative model is proposed and calibrated to available experimental data for auxetic LCEs. To analyse the effect of different model parameters on the predicted mechanical response, a simplified qualitative model is developed in [7]. The theoretical analysis suggests that auxeticity can be obtained by letting the uniaxial scalar order parameter decrease to a sufficiently low value then increase again. This is in agreement with the experimental observations where, in addition, mechanical biaxiality can also be observed that enhances the auxetic behaviour.
The purpose of this paper is to present a mathematical model for nematic LCEs that is capable of predicting a similar auxetic behaviour as that found in cat skin where reorientation of dermal fibres produces an increase in the thickness direction under tensile stretch. Cat skin thickening occurs first at small strain and continues at large strain until a maximum stretch is attained, after which the skin thickness begins to decrease. Both compressible and incompressible isotropic hyperelastic models are found acceptable within the range of experimental data, as reported in [15]. For the incompressible LCE model mimicking such behaviour, the uniaxial order parameter is allowed to decrease then increase again, and the critical stretch at which this change in monotonicity and where the director also rotates suddenly is predicted. In addition, a stochastic elasticity framework [5] is adopted where the model parameters are described by probability density functions, and the uncertainty in those parameters is propagated numerically to the predicted mechanical results. The theoretical model is described and explained in Section 2. This model is then calibrated to experimental data showing the auxeticity of cat skin in Section 3. The stochastic framework is summarised in the appendix.

Stochastic strain-energy function
In this paper, the model functions describing an incompressible nematic LCE satisfy the following assumptions [5]: (C1) Objectivity (frame-indifference). This condition states that the scalar strain-energy function is unaffected by a superimposed rigid-body transformation, which involves a change of position after deformation. As the nematic director n is defined with respect to the deformed configuration, it transforms when this configuration is rotated, whereas the nematic director n 0 for the reference configuration does not. (C2) Isotropy. This property requires that the strain-energy function is unaffected by a rigid-body transformation prior to deformation. As n is defined with respect to the deformed configuration, it does not change when the reference configuration is rotated, whereas n 0 does.
These fundamental assumptions underpin also the finite elasticity and liquid crystal theories [2,14]. In particular, the following two-term strain-energy function is considered, where μ 1 , μ 2 > 0 and m, n are given parameters, with μ = μ 1 + μ 2 representing the shear modulus at infinitesimal strain, {λ 2 1 , λ 2 2 , λ 2 3 } denote the eigenvalues of the Cauchy-Green tensor FF T , where F is the deformation gradient from the reference cross-linking state, satisfying det F = 1, and {α 2 1 , α 2 2 , α 2 3 } are the eigenvalues of AA T , where A = G −1 FG 0 is the local elastic deformation tensor, such that det A = 1, with G 0 and G representing the 'natural' deformation tensors due to the liquid crystal director in the reference and current configuration, respectively.
Given an intrinsically uniaxial LCE, the natural gradient tensor takes the form [7], where n is the nematic director, ⊗ denotes the tensor product of two vectors, I = diag(1, 1, 1) is the second order identity tensor, and represents the anisotropy parameter, with Q the uniaxial scalar order parameter. In the reference configuration, G is replaced by G 0 , with n 0 instead of n, a 0 instead of a, and Q 0 instead of Q. The model function defined by equation (1) can be regarded as a slight generalisation of the strain-energy density proposed in [7] where m = −n. Similar but more general models were introduced and analysed in [6]. In addition, here, it is assumed that the shear modulus μ follows a Gamma probability density function (see Appendix A for details), while m, n are deterministic constants.
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 LCE 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 θ ∈ [0, π/2] is the angle between n and n 0 , the deformation gradient takes the form Assuming 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 [6,7], we denote by λ crt the longitudinal stretch ratio λ 1 where the rotation occurs. In practice, although the angle θ may take intermediate values between 0 and π/2, the deformation interval for which the director rotation is observed can be very short, and its separate analysis will be omitted.
When the director for the reference and current configuration are given by (4), the associated natural deformation tensors are, respectively, The elastic deformation tensor A = G −1 FG 0 is then equal to In particular: • If θ = π/2, then For the corresponding elastic Cauchy-Green tensor A T A = diag α 2 1 , α 2 2 , α 2 3 , where "T " denotes the transpose (see also [6]), we have: • If θ = π/2, then 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 equal to Taking λ 1 = λ, we solve for λ 2 the equation To capture the auxetic response while maintaining the 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 ∈ (0, 1), at λ = 1, to a minimum value Q crt ∈ (0, Q 0 ) that occurs at the critical stretch ratio λ = λ crt , then increase again, also linearly, to a value Q 1 ∈ (Q crt , Q 0 ), i.e., such that Q(1) = Q 0 , Q(λ crt ) = Q crt and Q(λ max ) = Q 1 . Hence, When a is given by (3) and Q = Q(λ) is described by (17) , the energy function is equal to where W (lce) (λ 1 , λ 2 , λ 3 , a, θ) = W (lce) is defined by equation (1). Solving the equation gives the critical stretch ratio λ = λ crt .
In Fig. 1(a), the auxetic response, showing how the stretch ratio corresponding to the thickness direction first increases then decreases, is captured. In the plots, as λ 3 records a change of monotonicity at λ 1 = λ crt , a slight change in the gradient of λ 2 can also be observed, although this stretch ratio continues to decrease as λ 1 increases. Figure 1(b) depicts the associated deformation dependent uniaxial parameter Q given by equation (17). The energy function described by equation (18) is illustrated in Fig. 2(a). The corresponding first Piola-Kirchhoff stress defined by equation (15) is plotted in Fig. 2(b). The computer simulations were run in Matlab 2022a where inbuilt functions for random number generation were used.

Conclusion
Biomimetic materials are continuously sought in research fields such as biomedical engineering and regenerative medicine [1,3]. This paper presents a hypothetical LCE model that mimics the unusual elastic response of cat skin examined experimentally and modelled within the theoretical framework of isotropic finite elasticity in [15]. Specifically, when a material sample is extended longitudinally, assuming that its volume is preserved, its thickness increases until a critical deformation is reached, then decreases. The hypothetical model suggests that LCEs could also be produced that exhibit such auxetic response, both under the small and large strain regimes. The present study may also reignite interest in the extraordinary mechanical behaviour of cat skin, which has been less investigated and discussed than other skin tissue types [4].

Appendix A: Stochastic model parameters
For the stochastic modelling of nematic LCEs adopted in this paper, the following hypotheses are required [5,Chapter 6]: For any given finite deformation, at any point in the material, the shear modulus μ > 0 and its inverse, 1/μ, are second order random variables, i.e., they have finite mean value and finite variance. This assumption is guaranteed by setting the following mathematical expectations: The first constraint in (A.1) specifies the mean value for the random shear modulus μ, and the second constraint provides a condition from which it follows that 1/μ is a second order random variable. By the maximum entropy principle, the shear modulus μ with mean value μ and standard deviation μ = √ Var[μ] (which is equal to the square root of the variance, Var[μ]) follows a Gamma probability distribution with shape and scale parameters ρ 1 > 0 and ρ 2 > 0, respectively, such that The resulting probability density function takes the form where : R * + → R is the complete Gamma function The word 'hyperparameters' is used for ρ 1 and ρ 2 to distinguish them from μ and other material constants. When μ = μ 1 + μ 2 , assuming μ i > 0, i = 1, 2, we define the auxiliary random variable such that 0 < R 1 < 1. Then, the random model parameters can be expressed equivalently as follows, It is reasonable to assume in which case, the random variable R 1 follows a standard Beta distribution, with hyperparameters ξ 1 > 0 and ξ 2 > 0 satisfying where R 1 is the mean value and Var[R 1 ] is the variance of R 1 . The associated probability density function is for r ∈ (0, 1) and ξ 1 , ξ 2 > 0, (A.9) where B : R * + × R * + → R is the Beta function For the random coefficients given by (A.6), the corresponding mean values are and the variances and covariance take the form, respectively, Data Availability All data generated or analysed during this study are included in this article.

Conflict of Interests
The author has no competing interests to declare that are relevant to the content of this article.
Open Access This article is licensed under a Creative Commons Attribution 4.0 International License, which permits use, sharing, adaptation, distribution and reproduction in any medium or format, as long as you give appropriate credit to the original author(s) and the source, provide a link to the Creative Commons licence, and indicate if changes were made. The images or other third party material in this article are included in the article's Creative Commons licence, unless indicated otherwise in a credit line to the material. If material is not included in the article's Creative Commons licence and your intended use is not permitted by statutory regulation or exceeds the permitted use, you will need to obtain permission directly from the copyright holder. To view a copy of this licence, visit http://creativecommons.org/licenses/by/4.0/.